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Chapter  1 
Overview 


1.1  What  is  LSH 

LSH  is  a  code  for  analyzing  the  linear  stability  of  basic  states.  The  main  code  is 
generic  and  void  of  any  particular  physical  problem  yet  provides  for  various  tasks 
that  frequently  occur  in  the  stability  analysis  of  fluid  motions; 

global  —  eigenvalue  spectra, 

local  —  single  eigenvalues  and  eigenfunctions, 

table  —  one-  or  more-dimensional  tables  of  eigenvalues, 

curve  —  (neutral)  curves  in  parameter  space, 

and  others.  A  specific  physical  problem  can  be  defined  in  a  set  of  files  that  describe 
the  basic  state  and  the  stability  equations.  These  files  are  included  at  compile  time, 
while  definitions,  tasks,  and  peirameters  axe  read  from  two  additional  files  during  run 
time. 

The  original  version  C-S  of  LISA^  was  developed  to  be  able  to  obtain  some  addi¬ 
tional  results  on  some  stability  problem  in  a  few  minutes  if  the  insert  files  specific  to 
this  problem  are  available,  and  to  obtaiin  the  first  results  for  a  new  problem  within 
a  few  hours  if  the  equations  are  available  (see  LISA  Version  C-S  2.0  User’s  Manual). 
Since  the  specifications  for  a  new  problem  eire  usually  short,  the  debugging  effort  is 
greatly  reduced.  Often,  the  specification  of  a  new  case  can  be  eased  by  referring  to 
the  basic  state  or  the  stability  equations  of  an  existing  case  through  the  Makefile. 
This  Version  C-S  is  based  on  a  spectral  collocation  method  that  has  proven  reliable 
and  able  to  provide  highly  accurate  results  for  a  wide  spectrum  of  applications. 

The  present  Version  WL  of  LSH  is  a  derivative  of  LISA  and  has  a  similar  modular 
structure,  but  is  tailored  for  compressible  flows  over  aerodynamic  bodies  with  cur¬ 
vature.  Unlike  the  version  C-S,  the  stability  equations  here  axe  not  specified  by  the 
user.  They  are  already  coded  for  three-dimensional  disturbances,  derived  from  the 

^LISA  is  a  proprietary  software  product  of  DynaFlow,  Inc.  for  the  commercial  market. 
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Navier-Stokes  equations  for  compressible  flows.  The  analysis  is  done  in  general  curvi¬ 
linear  coordinates  defined  by  the  user.  The  basic  flow  to  be  studied  is  supplied  by 
the  user  either  by  analytical  calculation  or  reading  in  from  a  data  set.  The  user  may 
also  specify  boundary  conditions  for  disturbances  and  choose  thermodynamic  proper¬ 
ties  and  dependent  variables  for  disturbcinces.  Hence,  the  code  is  highly  modulair  and 
adaptable  to  different  problems.  The  stability  equations  in  this  version  are  discretized 
by  a  high-order  hermitijin  finite-difference  method.  This  method  is  accurate  yet  more 
efficient  than  the  spectral  method. 

1.2  What  is  PSH 

PSH  is  the  code  for  solving  the  compressible  PSE  for  three-dimensional  linear  or  non¬ 
linear  disturbances.  The  compressible  PSE  including  nonlinear  terms  are  formulated 
in  general  curvilinear  coordinates  and  are  already  coded  (Stuckert  &:  Herbert,  1993). 
The  user  should  refer  to  the  final  report  “Methods  for  Transition  Prediction  in  High¬ 
speed  Fows”  for  details  of  the  equations.  The  rest  of  the  physics  of  the  problem  is 
specified  by  the  user  by  insert  files  in  the  same  manner  as  for  LSH.  In  fact,  most  of 
the  generic  program  modules,  insert  files  and  an  input  file  Definitions  are  common 
to  PSH  and  LSH.  These  physics  files  are  described  in  detail  in  Section  3.  LSH  solves  an 
eigenvalue  problem  or  a  series  of  eigenvalue  problems.  Initial  conditions  for  PSH  are 
obtained  from  the  local  solution(s)  obtained  with  LSH.  The  PSH  then  integrates  the 
PSE  in  the  parabolic  marching  direction,  solving  the  system  (in  general,  nonlinear) 
of  equations  at  each  marching  step.  With  PSH,  linear  and  nonlinear  developments 
of  disturbances,  interaction  between  different  modes  of  disturbance  can  be  followed 
close  to  the  beginning  point  of  transition.  This  location  is  usually  meinifested  by  the 
minimum  in  mean  skin  friction  and  its  rapid  rise  thereafter,  etc.  Hence,  given  an  ac¬ 
curate  model  of  initial  data,  transition  can  be  predicted  without  any  empirical  data. 
The  computing  effort  required  for  PSH  is  much  less  than  that  for  Direct  Navier-Stokes 
simulations  (DNS). 


1.3  The  Program 

The  codes  LSH  and  PSH  are  highly  modulax.  They  are  structured  as  shown  in  Figure  1. 
For  this  Version  WL  1.0,  the  graphical  user  interface  (GUI)  is  under  development  for 
Silicon  Graphics  workstations  under  IRIX  4.0.5  or  higher.  Other  platforms  will  follow. 
The  batch  versions  have  been  tested  in  various  Unix  environments  and  can  be  adapted 
to  other  systems  by  minor  modifications,  especially  of  the  Makefiles.  The  codes  are 
partly  written  in  C  to  provide  the  connection  to  the  graphical  user  interface  and  for 
more  efficient  memory  management,  input,  and  output.  Other  parts,  including  the 
interface  to  the  user,  axe  written  in  Fortran  77  and  use  the  C  preprocessor  during 
object-code  generation.  The  following  describes  only  the  batch  versions  of  LSH  and 
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PSH. 

The  main  codes  are  provided  as  two  object  files,  Ish.o  and  psh.o,  eind  a  li¬ 
brary  archive  libut.a.  The  library  contains  the  eigenvalue  solvers  RG  and  CG  of 
EISPACK,  some  routines  of  LINPACK,  some  other  linear  algebra  solvers,  and  the 
routines  associated  with  the  spectral  collocation  method.  The  linear  algebra  routines 
consume  most  of  the  CPU  time.  All  floating-point  operations  are  performed  in  double 
precision,  real*8  or  complex*  16. 

In  spite  of  some  lacking  features,  Fortran  77  is  widely  available  and  used  in  en¬ 
gineering  applications.  The  goal  to  keep  the  code  generic  and  to  write  the  problem- 
oriented  “physics  insert  files”  in  Fortran  77  required  compromises.  Especicdly  the 
Fortran  hierarchy  of  statements  (which  prohibits,  e.g.,  dimension  after  executable 
statements)  made  the  placement  and  structure  of  the  “physics  insert  files”  manda¬ 
tory. 


Figure  1.1:  Prograun  Structure  for  LSH  and  PSH 
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Chapter  2 

The  Physics  Files 


Various  data  must  be  specified  to  adapt  the  generic  code  to  a  specific  problem.  This 
specification  is  made  through  two  groups  of  files:  insert  files  and  input  files.  There 
are  eight  insert  files  with  suffix  .ins  used  by  both  LSH  and  PSH,  i.e., 

common . ins 
gridbs . ins 
metric. ins 
bstate.ins 
btrans . ins 
boundary .  ins 
pointd.ins 


and  two  insert  files  with  suffix  .  ins  which  are  used  by  PSH  for  postprocessing: 

postpr . ins 
plotps . ins 

As  for  input  files,  there  is  one  common  to  both  LSH  and  PSH: 

Definitions 

one  input  file  specific  for  LSH: 

Parameters 

and  two  input  files  specific  for  PSH: 

Control 

Modes 
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The  insert  files  2ire  included  into  the  executable  code.  The  input  files  are  read  during 
run  time.  The  input  file  Modes  which  provide  initial  station  data  for  PSH  is  created 
by  running  LSH  at  the  desired  station  with  proper  parameters.  For  the  applications 
that  come  with  the  code,  the  user  has  to  deal  only  with  the  problem-specific  input 
files.  The  insert  files  axe  already  prepared.  For  new  applications,  the  developer  has 
to  prepare  the  insert  files  and  the  input  files.  This  preparation  usually  consists  of 
coding  the  calculation  of  the  b«isic  state  in  bstate.ins  and  filling  the  proper  array- 
elements  in  metric .  ins,  btrans .  ins,  etc. 


2.1  Directory  Structure  and  Msikefiles 

Before  we  discuss  preparation  and  contents  of  the  physics  files,  it  is  useful  to  consider 
the  arrangement  of  files  in  directories  and  how  the  executable  codes  for  a  specific 
stability  are  built.  Under  the  main  directory  w.psh  there  are  numerous  directories. 
The  first  character  of  the  directory  name  characterizes  the  contents: 

w.  *  :  general  objects  for  LSH  and  PSH 
c.* :  complex  problems 

b.*  :  basic  states;  gridbs .  ins ,  bst ate .  ins ,  btrains .  ins 
e.*  :  stability  equations;  common. ins,  metric. ins,  etc. 

All  specific  problems  are  in  separate  directories  c.  *.  Each  of  these  directories  contains 
at  leaist  four  files: 

Makefile 

Definitions 

Parameters 

Control 

Initially,  there  is  no  executable  code  Ish  or  psh  in  this  directory.  The  code  has  to  be 
made  using  the  Unix  commaind 


$  make 


By  default,  this  command  executes  the  instructions  in  the  Makefile.  The  first  line 
of  the  Madref  ile  defines  two  “macros”: 


ESTATE 

VECTOR 


=  directory  which  contains  gridbs. ins,  bstate.ins, 
and  btrans. ins; 

=  directory  which  contains  common,  ins  ,  metric,  ins, 
boundary .  ins,  and  pointd .  ins. 
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If  either  set  of  files,  say  bstate .  ins,  is  located  in  the  present  working  directory  c.case, 
then  either  BSTATE=... /c.case  or  BSTATE= would  be  appropriate  directives.  In 
the  latter  form,  the  left  quotes  request  that  the  system  set  BSTATE  to  the  output  of 
the  pwd  command.  This  latter  form  is  more  generic,  because  it  is  independent  of  the 
current  directory. 


Otherwise,  the  Makefile  leaves  the  task  of  making  the  executable  code  to  the 
.  ./Makefile  in  the  parent  directory  w.psh.  Figure  2.1  illustrates  how  the  Makefiles 
interact  with  the  directories  associated  with  the  problem  in  c.case.  In  essence,  the 
physics  insert  files  are  read  from  their  directories,  included  into  the  fortran  * .  f  files 
in  w.objs,  and  stored  in  w.tmp.  The  completed  files  in  w.tmp  are  compiled  into  object 
files  *.o  in  w.objs  and  combined  with  Ish.o,  psh.o  and  the  library  libut.a  to 
produce  Ish  and  psh  in  the  current  working  directory.  The  executable  codes  Ish 
and  psh  will  remain  there  even  if  the  user  works  on  other  problems,  until  the  file  is 
removed  with  the  Unix  command  nn  or 
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$  maike  clezin 


which  zJso  cleans  the  problem-specific  object  files  in  w.objs  and  prompts  for  the 
removal  of  the  output  files. 

Since  make  remakes  objects  only  if  they  are  out  of  date,  the  file  01d_dir  in  w.objs 
serves  to  register  if  the  user  changes  to  another  problem  with  different  insert  files. 
The  executable  code  Ish  or  psh  is  remade  if  the  insert  files  chainge. 

The  flexibility  in  accommodating  the  insert  files  in  the  same  or  different  directories 
avoids  unnecessary  duplications.  Often  it  is  only  convenient  to  protect  the  insert  files 
from  unintentional  damage  in  a  heavily  used  directory.  New  applications  must  be 
prepared  in  directories  on  the  same  level,  i.e.,  in  directories  w.psh/c.new-appl  to  be 
compatible  with  the  various  Makefiles. 


2.2  Preparing  the  Insert  Files 

In  the  following,  we  discuss  the  details  of  the  insert  files.  This  section  is  primarily 
provided  for  developers  of  new  applications.  All  insert  files  have  access  to  the  following 
shared  information: 

zro  =  0,  the  real  double  precision  number  zero 

one  =  1,  the  real  double  precision  number  one 

zero  =  (0,0),  the  complex  double  precision  number  zero 

uimag  =  (0,1),  the  complex  double  precision  imaginary  unit 

pi  the  real  double  precision  number  r  =  3.1415  •  •  •. 

jx  the  highest  index  of  the  collocation  points,  j  =  0,  •  •  •  ,jx. 

nx  the  highest  order  of  differentiation. 

neon  the  number  of  constants  or  zero  (if  there  are  no  constants), 
npar  the  number  of  parameters. 

nequ  the  number  of  equations  in  the  stability  equations, 

nvar  the  number  of  variables  nvar  =  nequ  in  the  problem, 

isize  the  dimension  {jx  -I- 1)  •  nequ  of  the  algebraic  system. 
con(i)  the  constants  Ci  for  i  =  1,  •  •  • ,  neon. 

parCi)  the  parameters  Pjt  for  A:  =  1,  •  •  • ,  npar. 

chy  ( j  )  the  collocation  points  (j,  j  =  0,  •  •  •  ,jx. 

The  following  variables  related  to  PSH  are  also  available. 

modt  the  highest  index  of  the  temporal  frequency,  the  lowest  index 

is  0. 

modzm  the  lowest  index  of  wavenumber  in  direction 
modzp  the  highest  index  of  wavenumber  in  direction 
izsym  integer  indicating  symmetry  direction  for  modes 


7 


pann(k,it,iz,l)  the  parameters  Pk  for  k  =  l,...,npar;  and 

modes  it  =  0,  ...,modt-,  iz  —  —  izsym)  * 

modzm, modzp  at  /  =  Iprev  or  /  =  Icur  step. 
yofx(j  ,1)  the  collocation  points  j  =  0,  •  ■  •  at  /  =  Icur  or  I  = 

Iprev  step. 

dydx(j )  at  j  =  0,  •  •  • ,  ji  at  the  current  step. 

The  variable  names  follow  the  Fortran  convention  that  all  names  beginning  with 
i-n  are  implicitly  integer.  All  variables  with  names  beginning  with  a-h  or  o-z  are 
implicitly  declared  double  precision,  real*8.  Most  of  these  variables  tire  taken  from 
the  input  files  Definitions,  Parameters,  or  Control. 


2.2.1  common.ins 

This  file  is  for  communication  between  insert  files  bstate .  ins,  btrans .  ins,  metric .  ins, 
etc.,  and  the  stability  equations  of  the  main  code.  The  file  allocates  local  arrays  for 
all  basic  state  quantities,  metric  terms  and  declares  all  complex  variables,  used  in 
the  stability  equations.  The  maximum  number  of  basic-state  collocation  points  per 
station  to  be  read  in,  Fortran  I/O  unit  numbers  for  input  files  associated  with  basic 
state  and  user- written  output  files  are  also  specified  with  parameter  statements.  The 
specification  of  jmax  is  independent  of  jx  used  for  stability  analysis. 

No  executable  statements  are  allowed  in  common.ins. 

2.2.2  gridbs.ins 

Subroutine  gridbs  generates  the  grid  for  the  b«isic  state  calculation.  It  usually  does 
this  simply  by  reading  in  the  grid  of  a  previously  computed  basic  state. 

For  instance,  if  the  basic  state  is  computed  with  a  separate  boundary  layer  code, 
PNS  code,  or  Navier-Stokes  code,  then  gridbs  must  read  in  the  corresponding  grid. 
This  grid  should  then  be  made  available  through  a  common  block  to  subroutine 
bstate  which  will  have  to  interpolate  the  solution  onto  the  disturbance  grid. 

If,  however,  the  basic  state  solution  is  analytical  or  can  be  otherwise  easily  com¬ 
puted,  the  chy(j)  should  be  used  as  the  basic  state  grid.  In  this  case,  gridbs  is 
redundant,  and  the  file  gridbs.ins  may  be  empty  since  chy(j)  is  directly  available 
to  bstate. 

Subroutine  gridbs  may  also  rescale  the  basic  state  grid  and  solution  with  the  dis¬ 
turbance  scales,  or  vice  versa,  and  perform  other  necessary  basic  state  initializations. 

Like  bstate. ins  in  bstate. f,  gridbs.ins  is  included  in  gridbs. f  before  any 
excutable  statement.  Hence,  the  user  may  declau-e  new  arrays  as  needed. 
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2.2.3  metric.ins 


The  coordinate  system  chosen  for  the  stability  ajialysis  must  be  specified  by  the  user 
by  2issigning  values  to  the  covauiant  basis  vectors  dx‘/d^^  and  their  first  and  second 
partial  derivatives  with  respect  to  in  the  file  metric .  ins.  Here  x'  is  the 

Cartesian  coordinate  and  is  the  curvilinear  coordinate  used  in  the  stability 

analysis.  Hence,  the  following  corresponding  airrays  axe  defined  in  this  insert  file: 

dxdc{2 1  i)  =  dx* / ,  covaxiant  basis  vectors  =  l,nd  where 
nd  is  the  number  of  dimensions  .  nd  is  4  usually,  with  3 
spatial  coordinates  and  a  fourth  temporal  coordinate. 

d2xdc2(k,  j ,  i)  =  first  partial  derivatives  of  covariant 

basis  vectors  with  respect  to  k  =  l,nd. 

d3xdc3(l,k,  j, i)  =  second  partial  derivatives  of 

covariant  basis  vectors  with  respect  to  and  fc  =  l,nd;, 

/  =  l,nd. 

Symmetry  of  metric  terms  axe  used  in  the  subroutine  metric  so  that  the  user  needs 
to  specify  only  unique  entries  of  arrays  d2xdc2  and  d3xdc3.  The  array  d3xdc3  needs 
not  be  specified  if  Cartesian  velocity  components  are  used  as  dependent  variables  for 
stability  analysis  (mvar=0  is  specified).  Subroutine  metric  uses  these  in  metric.ins 
to  compute  the  terms  d^'^jdx^,  jdx^dx^,  and  fdx^dx’^dx^  which  appear  in 
the  stability  equations  after  transforming  the  x‘  derivatives  of  the  disturbances  to 
derivatives  with  the  chain  rule.  It  is  important  to  note  that  the  approximations  to  the 
stability  equations  (either  local  or  PSE)  axe  made  after  the  coordinate  transformation. 
Hence,  the  coordinate  system  chosen  will  have  some  effect  on  the  accuracy  of  the 
results  obtained. 

For  each  streamwise  station  ,  these  covcuriant  and  contravariant  metric  terms  for 
all  {jx)  normal  grid  points  are  saved  in  global  arrays.  These  arrays  are  used  in  other 
modules  of  the  code  when  necessary.  For  every  n  —  th  collocation  point  at  a  station, 
the  following  metric  terms  are  available  for  the  user,  e.g.,  to  postprocess  the  station 
solution  in  the  insert  files  postpr.ins. 

dxdca(n, j ,i) 
d2xdca(n,ks,i) 
d3xdca(n,ls,i) 


gda(n) 


dx' 

w 

a^x* 

dx^  dx^ 

ww 


9 


sqrtga(n) 


dcdxa(n,j ,i) 

d2cdxa(ii,ks,i) 

d3cdxa(n,ls,i) 


1 

\  d(' 

Jacobian 

dxi 

dH' 

dxidx’^ 

dH' 

dx^dx^dx^ 


t  =  1,  ...,nd 

{{{{Is  =  /5  +  /),  /=  ks  =  ks  +  k),k  =  1, 


2.2.4  btrans.ins 

In  stability  analysis,  the  equations  are  derived  for  disturbances,  Q',  i  =  l,...,nvar 
about  a  laminar  basic  flow,  Q*.  The  stability  code  requires  that  basic  state,  in  vari¬ 
ables  Q'  be  specified  for  each  grid  point  in  the  stability  analysis  coordinate  system 
specified  in  metric. ins.  The  basic  state  is  usually  obtained  in  the  insert  file 
gridbs.ins  bstate.ins  in  variables  and  a  coordinate  system  suitable  for  its  solu¬ 
tion.  Instead  of  requiring  the  user  to  tramsform  these  into  Q'  in  stability  coordinates 
directly,  more  flexibility  is  provided  by  using  the  subroutine  btrans.f  which  allows 
specification  of  three  sets  of  variables  through  insert  file  btrans .  ins.  Note  that  sub¬ 
routine  btrans  heis  executable  statements  before  including  btrans .  ins.  Therefore, 
no  array  declarations  should  be  made  in  this  insert  file.  The  insert  file  btrains .  ins 
must  be  coded  so  that  it  performs  the  following  tasks. 

1.  The  variables  defined  in  bstate.ins  are  transformed  into  density  and  temper¬ 
ature  as  the  first  nts  =  2  dependent  variables  and  one  of  the  following  sets  of 
velocity  components  as  nts  -I- 1  to  nvar  variables: 

(a)  Cartesian  velocity  components. 

(b)  Physical  contravariant  velocity  components  in  the  stability  coordinate. 

(c)  Contravariant  velocity  components  in  the  stability  coordinate. 

The  user  must  specify  one  of  the  above  three  sets  by  setting  ibvel  to  either 
0, 1,  or  2,  respectively. 


2.  The  derivatives  of  the  variables  obtained  in  bstate.ins  are  transformed  into 
derivatives  with  respect  to  curvilinear  stability  coordinates  . 

3.  The  derivatives  of  the  equation  of  state  and  specific  enthalpy  with  respect  to 
thermodynaunic  dependent  variables  are  assigned. 

Hence,  the  following  arrays  of  basic  state  variables  are  to  be  defined  by  the  user 
in  this  insert  file. 

bsoln(i)  =  Q\  i-th  dependent  variable  for  basic  state,  with  2  =  1 
to  nts  representing  thermodynamic  dependent  variables  p  and 
T,  and  nts  +  1  to  nvar  Cartesian,  physical  contravaxiant  or 
contravariant  velocity  components. 

dbdcCj  ,i)  =  dQ'fd^^,  first  partial  derivatives  of  Q*  with  respect 
to  curvilinear  coordinates,  =  l,nd  where  nd  is  the  num¬ 
ber  of  spatial  dimensions,  usually  3. 

d2bdc2(k,  j  ,i)  =  d'^Q' ,  second  partial  derivatives  of  basic 
state  with  respect  to  and  j  =  l,nd;  k  =  l,nd. 

ddqlp(i),  ddq2p(i,  j),  ddq3p(ig,k),  first,  second,  and  third  par¬ 
tial  derivatives  of  the  equation  of  state  p{p,  T)  with  respect 
to  thermodynamic  variables,  Q\  i  =  I,  nts,  i.e., 

^P  •  •  i.  _  1  t 

dQi'  dQWQ'^  dQ^dQ^dQ''  ' 

For  example,  for  a  perfect  gas,  p  =  Hence 


ddqlp(l) 

ddqlp(2) 

ddq2p(l,l) 

ddq2p(l,2) 

ddq2p(2,2) 


_  ^  _  T 
~  dp 

_  ^  _  P 
dT 
=  0 

=  ddq2p(2,l)=^ 
=  0  etc. 


ddqlh(i),  ddq2h.(i,  j),  ddq3h(i,j  ,k),  first,  second,  and  third  par¬ 
tial  derivatives  of  specific  enthalpy  h{p,  T)  with  respect  to 
thermodynamic  variables  ,  i.e., 

dk  d^k  m  •  •  ,  1 

ag*’  ag^agw’ 
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ddqlml(i),  ddq2ml(i,j),  first  and  second  partial  derivatives  of 
first  coefficient  of  viscosity  ni  with  respect  to  thermodynamic 
variables  ,  i.e., 


dQi'  dQ^dQ''  ~ 

ddqlm2(i),  ddq2m2(i, j),  first  and  second  partial  derivatives  of 
second  coefficient  of  viscosity  ^2  with  respect  to  thermodynamic 
variables,  i.e.  , 


dQ''  dQidQ^'  ~ 

ddqlk(i) ,  ddq2k(i ,  3 ) ,  first  and  second  partial  derivatives  of  ther¬ 
mal  conductivity  «  with  respect  to  thermodynamic  variables, 
i.e.. 


Bk  d^K  .  . 
dQidQ''  = 


Higher  order  partial  derivatives  of  thermodynamic  coefficients  with 
respect  to  temperature.  These  are  needed  for  evaluating  non¬ 
linear  terms. 


dSmuldt 

d4iauldt 

d3kdt 

d4kdt 

d4hdt 

dBhdt 


d*tii 

dT* 

d^K 

d*K 

dT* 

d*h 

dT* 

d^h 

dT^ 


Note  that  the  names  of  variables  here  do  not  necessarily  coincide  with  those  used 
in  bstate.ins  for  the  basic  flow  solution.  Assignment  statements  in  btrans.ins 
should  translate  the  names  of  the  latter  into  those  of  the  former.  For  exaimple, 
bstate.ins  might  compute  the  similarity  solution  for  the  compressible  boundary 
layer  and  have  as  output  the  streamfunction  and  its  derivatives  with  respect  to  simi¬ 
larity  coordinates.  The  user  would,  thus,  transform  the  similarity  variables  into  one 
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set  of  velocity  components  (described  above)  in  the  file  btrans .  ins  axid  transform 
derivatives  with  respect  to  similarity  coordinates  into  those  with  respect  to  curvilinear 
stability  coordinates.  ^ 

Using  these  user  specified  basic  state  values  and  metric  terms,  the  subroutine 
btrans  calculates  the  following  quantities  in  Cartesian  coordinates  and  stores  them 
in  global  arrays.  They  are  available  at  all  n  collocation  points  for  use  in  the  insert 
files. 


bsolna(n,i) 
dbdxa(n, j ,i) 

d2bdxa(n,ls,i) 

blapa(n,i) 

bgrddv(n,j) 


Q‘  =  {p,T,u,v,wy 

dQ' 

dxi 

dx^dx^ 

AQ 

d  du'dx' 
dxi 


The  following  thermodynamic  quantities  ar 

thermka(n)  = 

amul (n)  = 

ainu2  (n)  = 

dlqk(n,i)  = 

d2qk(n,ls)  = 

dlqmuKn,!)  = 

d2qmul (n,ls)  = 

dlqmu2(n,i)  = 

d2qinu2(n,ls)  = 

dlqh(n,i)  = 

d2qh(n,ls)  = 

^One  might  as  well  do  this  directly  in  bstate. 
the  basic  state  is  often  independent  of  the  basic  stat 
it  is  easier  to  modify  the  basic  state  numerics  whe 


also  available  at  all  n  collocation  points. 


I^i 

Bk 

W 

d^K 

BQ'BQi' 

d}h_ 

BQ'^ 

B^Pi 

BQ'BQC 

dp2 

oQ'' 

d^P2 

BQ'BQ^' 

Bh 

BQ'' 

B^h 

BQiBQ^' 

ns.  However,  this  preprocessing  (translation)  of 
solution  process.  With  this  separate  trainslation, 
necessary. 
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d3qh(n,ks) 


dlq-thermo  (n ,  i) 


d2q_thermo  (n ,  Is) 


d3q>thermo  (n  ,ks) 


d^h 

dQ'dQ^dQ>^' 

d{state  equation} 

’ 

d^{state  equation} 
dQ^dQ^  ’ 
d^{state  equation} 
dQ'dQ^dQ‘‘  ’ 


((((fcs  =  0,...,ks  +  k),  k  =  Is  =  0, 


Is  +  j),  j  =  1, 


i),  i  =  1,  ...,n<s) 


Note  that  second  and  higher  derivatives  axe  symmetric,  i.e., 


_f(i_  =  J!ii_ 

dQ'aQ^  dQ^dQ'’ 

Hence,  only  unique  entries  of  the  derivatives  are  stored  in  these  arrays. 


2.2.5  pointd.ins 

After  the  coordinate  system  has  been  specified,  the  grid  point  distribution  for  the 
stability  analysis  must  be  chosen.  This  is  done  partly  by  specifying  the  choice  of 
stretching  function  idtrf ,  and  parameters  (strfp(i)  ,  i=l , ,4)  in  the  input  file 
Definitions  and  partly  through  the  insert  file,  pointd.ins.  The  stability  grid  is 
given  cis  Different  choices  of  idtrf  described  in  Section  2.3 

specify  insert  file  pointd .  ins  is  included  in  the  subroutine  f  dinit  where 

the  stretching  function  is  evaluated.  In  pointd.ins,  the  user  may  specify  how  the 
grid  point  distribution  in  the  direction  changes  with  distance  downstream,  i.e. 

For  example,  the  location  of  the  outer  boundary,  (stored  in  the  local 
variable  btrf )  may  be  varied  as  a  function  of  the  streamwise  distance  to  prescribe; 


e{j,e)=eef 


The  value  of  at  the  reference  station  is  input  as  strfp(2).  If  pointd.ins  is 
empty,  then  the  same  grid  point  distribution  in  at  reference  location  will  be  used 
for  the  analysis. 

^  Here,  different  choices  do  not  affect  the  accuracy  of  the  approximations  made  to  the  stability 
equations,  but  do  affect  the  numerical  accuracy  with  which  the  discretized  equations  are  solved. 
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2.2.6  boundary.ins 

In  this  insert  file,  boundary  conditions  for  disturbances  are  specified.  This  insert  file 
is  included  in  the  subroutine  bvector.  The  user  can  choose  the  dependent  variables 
/  for  the  stability  ^alysis  by  using  the  flag  mvar  in  the  input  file  Definitions.  The 
default  variables  (Q,  mvax=0)  are  p,  T,  it*,  where  p  is  the  density  disturbance,  f  is  the 
temperature  disturbance,  and  u*  is  the  disturbance  in  the  i-th  Cartesian  component 
of  the  disturbance  velocity.  Transformation  from  /  to  Q  is  described  by  the  following 
equation. 

Q  =Mf  or  Q'  =Mir 

This  matrix  M  and  its  derivatives  with  respect  to  Cartesian  coordinates,  for  a  speci¬ 
fied  value  of  mvar,  are  evaluated  in  the  subroutine  setvrm.  They  are  available  in  the 
following  global  real  arrays  at  n  =  0, ...,  jx  collocation  points  for  the  user. 


vrma(n, j ,i) 
dlxvraCn.k, j ,i) 

d2xvra(n,ks, j ,i) 


M'j  ;  2  =  1, . . . ,  nvar^j  =  1, . . . ,  nvar 


dM] 

W  ’ 


A;  =  1,...,4 

;  (((A:s  =  0,ks  +  l),l  =  1,. 


..,3) 


Global  arrays  for  metric  terms  and  basic  state  (as  described  in  metric. ins  and 
btrans .  ins)  are  also  available  here. 

The  boundary  conditions  axe  applied  on  the  variables  /  chosen  for  stability  analy¬ 
sis.  In  general,  the  boundary  condition  for  the  f-th  equation  is  imposed  by  specifying 
the  real  array  bvec(i,m,l,k)  in  the  boundary-condition  equation 


dmfl 

bvec(i,m,l,k)^^ 


bvrhs(i,k)  ;  I  =  1, . . . ,  nuar,  2  =  1,...,  nequ 


Here,  m  =  0,  ...,nx.  A:  =  0  or  1  for  the  boundary  condition  at  or 

respectively.  An  integer  array  numbc(i,k)  should  be  set  to  the  number  of  boundary 
conditions  specified  (0  for  no  conditions). 

In  many  problems,  physical  quantities  such  as  Q^=p,  Q^=T  and  contravariant 
velocities  (not  /•'  ),  are  readily  known  at  the  boundaries.  Then  appropriate  trans¬ 
formations  may  be  used  in  the  boundary  conditions.  When  Neuman-type  boundary 
conditions  or  mixed  Neuman-Dirichelet  boundary  conditions  on  Q^-,  j  =  1,2  are 
known  at  the  boundaries,  derivative  of  Q  with  respect  to  may  be  related  to  that 
of  /  as  follows. 

dQ  dM 

de  “  de  ^ 
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dM can  be  evaluated  from  the  stored  global  2irrays  dlxvra(n  ,k ,  j ,  i)  and  dxdcaCn ,  2  ,k) 
by  using  the  chain  rule. 

dMj  dM]  dx‘‘ 

de  ~  W 


Similaxly,  boundary  conditions  on  contravariant  velocities  u^”*  may  be  related  to  those 
of  p  ais  follows: 


-C" 


.  ,3;  m  =  1, . . .  ,3 


_  For  example,  one  may  specify  the  adiabatic  boundary  condition  on  f  (t  =  2),  i.e., 
dT fd(^  =  0  at  the  wall  (A:  =  0)  as  follow: 


bvec(2, 1,1,0) 
bvec(2, 0,1, 0) 
bvrhs(2,  l) 


vrma(0, 1, 2) 
dlcvrm(2, 1, 2) 
0 


/  =  !,...  ,nvar 


n\mbc(2, 0)  =  1 

Here,  dlcvrm(2,l,2)  is  the  local  aurray  for  dMfldi"^  evaluated  at  the  collocation 
point  n  =  0. 

The  following  general  conditions  aire  applied  in  boundary .  ins  files  provided  with 
the  built-in  sample  applications.  At  the  wall,  disturbance  temperature  and  ve¬ 
locities  aire  set  to  zero.  At  farfield,  disturbance  temperature  velocities  in 
and  directions  are  set  to  zero.  Then  ^^-gradient  of  disturbance  velocity  in 
direction  is  determined  from  the  continuity  equation  assuming  that  density  and 
gradients  of  density  are  zero  implicitly.  ^  No  explicit  boundary  condition  on  den¬ 
sity  is  required  at  either  end  of  the  boundaries.  For  the  current  version  of  the  code, 
aisymptotic  conditions  (Appendix  B)  au:e  implemented  only  for  blunt-cone  calculations 
using  orthogonal  body-intrinsic  coordinates.  Boxmdary  conditions  for  the  stationary 
{u>  =  0)  modes  may  be  different  from  those  of  oscillatory  disturbamces.  For  example, 
when  the  basic  laminar  flow  satisfies  the  auliabatic  wall  boundary  condition,  normal 
gradient  of  T  for  stationary  modes  may  be  set  to  zero  (see  Mack  1984). 


2.2.7  postpr.ins  for  PSH 

This  insert  file  is  for  postprocessing  the  data  at  specified  stations  to  obtain  specific 
output  quantities.  This  file  is  included  in  subroutine  wmodes  before  any  executable 

^The  user  may  choose  to  set  this  velocity  component  to  zero  for  oscillatory  modes  with  little  or 
no  change  in  the  solution. 


16 


statement.  Hence  the  user  may  declare  new  variables  at  the  top  of  this  insert  file,  as 
required.  At  each  streamwise  station,  the  user  has  access  to  the  following  arrays: 

func(n,m,iv,it,iz,lt)  =  ir-th  variable  of  the  shape 

function  f  for  mode  {it,iz)  and  its  m-th  derivatives  with  re¬ 
spect  to  at  n-th  collocation  point  and  {It  =  Iprev)  previous 
step  or  {It  =  Icurr)  current  step  (complex), 
n  =  0, ...,  jx;  m  =  0,  ...,nx;  iv  =  1,  ...,nvar. 

exprCit  ,iz,lt)  =  exp~^'^\  amplitude  modulation  owing  to  ex¬ 
ponential  growth  for  mode  {it,iz)  at  {It  =  Iprev)  or  {It  = 

Icurr)  step,  (complex). 

paxm(k,it  ,iz,lt),  parameters  P*,  h  =  1,  ...,npar,  such  as  wavenum¬ 
bers  and  frequency  (as  defined  in  definitions)  for  mode 
{itjiz)  at  the  current  or  previous  step  (real). 

The  user  may  also  use  the  available  basic-state  cirrays  described  in  btrans .  ins, 
metric  arrays  described  in  metric .  ins,  and  variable  transformation  airrays  described 
in  boundary .  ins. 

2.2.8  plotps.ins  for  PSH 

This  insert  file  is  included  in  subroutine  wmodes  after  postpr.psh.  It  opens  two 
output  files  for  a  mode  {it,iz),  as  this  mode  is  introduced. 

1.  Profile  data  file  pt(it)zp(iz)  .dat 

2.  Station  data  file  st(it)zp(iz)  .dat, 

where  it  and  iz  in  parentheses  are  to  be  substituted  by  indices  of  the  mode.  The 
user  is  encouraged  to  use  the  same  logical  statements  in  naming  these  output  files, 
although  other  contents  of  these  files  also  can  be  changed. 

The  Profiles  of  eigenfimctions  and  growth  rates  eire  written  to  pt(it)zp(iz)  .dat 
file  at  every  npprt  step  (station).  The  value  of  npprt  is  specified  as  a  constant  in  the 
Definitions  file. 

Station  (current  step)  data  are  printed  out  to  st(it)zp(iz)  .dat  every  step  in 
It  may  include  parameters,  maximum  amplitude  of  different  disturbance  vari¬ 
ables,  growth  rates  along  constant  amplitudes  of  meissflow  and  total  temperature 
fluctuations  at  the  location  where  the  maissflow  for  the  reference  mode  {itr,izr)  is 
maximum,  and  their  corresponding  growth  rates,  etc. 

Some  of  the  station  data  are  also  printed  out  in  0 .  summary,  the  standard  output 
file  for  PSH.  Insert  files  postpr.psh  and  plot  .psh  are  provided  so  that  the  user  may 
customize  the  output  for  a  particular  problem. 


2.3  Preparing  the  Definitions 

The  file  Definitions  is  read  in  a  format  similar  to  the  free  format  in  Fortraui  77. 
Exceptions  are:  all  items  must  be  separated  by  spe  ’es  (no  commas  are  allowed):  the 
items  may  be  character  strings;  and  the  data  cannot  be  continued  on  the  next  line. 
The  slash  /  terminates  the  reading  of  data  from  the  line.  Empty  lines  are  ignored 
and  can  be  freely  used  to  enhance  readability. 

The  file  Definitions  consists  of  various  groups  of  data  described  in  the  following: 

(1)  A  single  line  of  text  specifying  the  problem. 

(2)  One  line  specifying  the  integer 

jx  the  highest  index  of  the  collocation  points,  j  =  0,  •  •  •  ,jx.  The 
number  jx  +  1  of  points  can  reach  from  as  few  as  30,  say,  to 
as  many  as  200  or  so. 

(3)  One  line  specifying  the  integer 

neon  the  number  of  constants  or  zero  (if  there  are  no  constants). 

(4)  One  line  specifying  the  integer 

the  number  of  parameters.  Note  that  complex  parameters 
count  as  two  real  parameters. 

(5)  One  line  specifying  the  type  of  collocation  points  for  stability  analysis,  idpts: 

1  for  gauss  points  in  (—1, 1)  for  spectral  version  only. 

0  for  using  the  b2isic-state  grid  points 
-1  for  stretched  grid  points 

Normally,  uniformly  spaced  grid  points  are  suitable  for  Hermitian  finite  dif¬ 
ference  method  used.  For  high-speed  boundary  layers,  the  grid  resolution  is 
required  not  only  near  the  wall  but  also  near  the  boundary  layer  edge  where 
the  critical  layer  is  located. 

(6)  One  line  specifying  the  stretching  in  (^{j)  with  the  following  five  data: 
idtrf  the  type  of  transformation  (integer) 

strfp(l)  the  minimum  value  of  a,  (real) 

strfp(2)  the  maximum  value  of  b,  (real) 

strfpO)  the  third  stretching  parameter,  or  zero  (real) 
strfp(4)  the  fourth  stretching  parameter,  or  zero  (real) 

The  following  transformations  are  avaulable  by  specifying  idtrf: 

1  Equally  spaced  grid  (^(j)  =  a-\-{b  —  a)*j fjx 
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2  Exponential  stretching  ^*(j)  =  a  -t-  (6  -  a)  *  (H  -  l)/(r-’^  -  1) 
where  strfp(3)=r(''®~^)  is  the  ratio  of  maximum  grid  spacing 
at  (^ax  minimum  spacing  at  strfpO)  must  be  greater 
than  or  equal  to  1. 

3  Exponential  stretching 

eu)  =  <■ + (t  -  g) » ^  - 1} 

^strfp{3)  - 1^ 

4  Logarithmic  stretching  of 

(^(j)  =  a  +  (6  -  a)  *  (ru/rl) 
with 

ru  =  log{strfp(Z) +j/jx)/{strfp{3)  -  j/jx) 
rl  =  log{sirfp{3)  +  l)/(str/p(3)  -  1) 


5 


Typically,  strfpO)  «  1.1. 
The  hyperbolic  sine  mapping 


eu) 

where 


a  +  {3tr/p(3)  -  a} 


I  +  sinh{strfp{i))j^ 
sinh{strfp{A)  c) 


-  strfp{4)log 


1  ^  _  l)»*’‘/P(3)-a 

1  -I-  (e-»t’-/p(4)  - 


strfp{3)  >  a,  strfp{A)  >  0 


(7)  nvar  lines,  where  each  line  specifies  one  variable  by  the  following  data: 
vname  the  ncune  of  the  disturbance  variable  p,  up  to  8  characters, 
iiegv  the  ^^-symmetry  of  this  variable,  where 

0  indicates  no  symmetry  condition  on  this  variable 
1  indicates  this  vjiriable  is  symmetric 

-1  indicates  this  variable  is  ^ultisymmetric 

idsy  the  ^^-symmetry  of  this  variable,  where 
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0  indicates  no  symmetry  condition  on  this  variable 
1  indicates  this  variable  is  symmetric 

-1  indicates  this  viiriable  is  antisymmetric 

(8)  neon  lines,  where  each  line  specifies  one  constant  by  the  following  data: 
cname  the  name  of  the  constant  Ci,  up  to  8  characters. 

con  the  value  of  the  constant  Ci  (reaJ) 

This  section  may  be  missing  if  ncon=0. 

(9)  npar  lines,  where  each  line  specifies  one  paxauneter  by  the  following  data: 
pname  the  name  of  the  parameter  Pk,  up  to  8  characters. 

ipar  the  appearance  of  the  parameter  in  the  problem,  where 

0  indicates  that  the  parameter  appears  linear  in  the  sta¬ 
bility  equations  and  does  not  affect  the  basic  state 

1  indicates  that  the  parameter  appears  nonlinear  in  the 
stability  equations  and  does  not  affect  the  basic  state 

2  indicates  that  the  parameter  affects  the  basic  state 
it3rp  the  real  or  complex  nature  of  the  parameter,  where 

0  indicates  that  the  parameter  is  real 

1  indicates  that  the  paraimeter  is  the  real  part  of  a  com¬ 
plex  quantity 

2  indicates  that  the  parameter  is  the  imaginary  part  of  a 
complex  quantity 

Real  amd  imaginary  parts  of  complex  quantities  must  immediately  follow  each 
other. 

Although  tedious  to  describe,  the  Definitions  are  usually  quick  done.  However,  it 
is  important  to  prepare  this  file  carefully  since  the  sequence  of  data  must  be  correct. 
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Chapter  3 

Tasks  and  Parameters  for  LSH 


The  file  Parameters  contains  the  sequence  of  tasks  for  the  code  LSH  to  perform. 
The  file  Parameters  is  read  in  a  format  similar  to  the  free  format  in  Fortran  77. 
Exceptions  are:  all  items  must  be  sep^ated  by  spaces  (no  commas  are  allowed);  the 
items  may  be  character  strings;  and  the  data  cannot  be  continued  on  the  next  line. 
The  slash  /  terminates  the  reading  of  data  from  the  line.  Empty  lines  are  ignored 
and  can  be  freely  used  to  enhance  readability. 

The  code  can  perform  the  following  tasks: 


it  ask  =  1 
it  ask  =  2 
itask  =  3 
itask  =  4 

itask  =  -1 
itask  =  -2 
itask  =-999 


local  procedure  to  find  a  single  eigenvalue 
local  procedure  to  find  eigenvalue  and  eigenfunction 
table  of  eigenvalues  in  1,  2,  or  3  dimensions 
curve  of  solutions  in  the  plane  of  two  parameters 

close  and  reopen  the  file  Parameters  after  pause 

execute  a  system  command 

terminate 


These  tasks  and  their  required  input  are  described  in  the  following  sections. 


3.1  Single  Eigenvalues 

Procedure  to  find  spectra  of  eigenvalues  is  not  yet  implemented  in  the  Hermitian  ver¬ 
sion  of  LSH.  In  contrast  to  the  global  procedure  to  find  spectra,  the  local  procedure 
searches  for  the  value  of  one  real  parameter  (for  real  problems)  or  the  values  of  two 
real  parameters  (for  complex  problems)  that  solve  the  characteristic  equation.  For 
complex  problems,  the  two  real  parameters  Me  not  necessarily  the  real  and  imaginary 
part  of  a  complex  quantity.  We  denote  these  two  parameters  as  the  “first  and  second 
eigenvalue  parameter.”  To  obtain  an  eigenvedue  with  the  local  procedure,  an  initial 
guess  must  be  specified.  If  such  an  estimate  is  not  available,  itask=0  can  help  to  pro¬ 
vide  many  eigenvalues  (if  the  eigenvalue  appears  linearly  in  the  stability  equations). 
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Otherwise,  itask»3  (table)  can  be  used  to  analytically  continue  an  eigenvalue  known 
at  different  parameters  to  the  desired  point. 

To  obtain  a  single  eigenvalue  for  a  certain  parameter  combination,  the  file  Parameters 
requires  the  following  information: 


(a)  One  line  with  it  ask: 

1  to  request  an  eigenvalue 


(b)  One  line  with  the  following  data: 

ievl  the  index  of  the  first  eigenvalue  parameter 

iev2  the  index  of  the  second  eigenvalue  parameter  (for  complex 

problems  only) 

-1  an  optional  specification  to  reset  the  eigenfunction 


(c) 


npar  lines,  where  each  line  specifies: 

par(k)  the  value  of  parameter  Pk.  The  values  of  par  (ievl)  and 
fortpar(iev2)  are  set  to  the  initial  guess. 


For  the  first  eigenvalue,  the  iteration  uses  the  initial  guess  for  the  eigenvzdue  and  a  de¬ 
fault  setting  for  the  eigenvector.  For  subsequent  eigenvalues,  the  converged  eigenvalue 
and  eigenvector  of  the  previous  solution  is  used  as  an  initial  guess  to  speed  up  the 
convergence.  The  optional  reset  of  the  eigenvector  to  the  default  setting  breaks  this 
routine.  This  option  is  provided  to  avoid  problems  with  changes  from  one  eigenmode 
to  another,  e.g.,  from  an  antisymmetric  to  a  symmetric  mode.  This  option  should  be 
used  if  the  mode  changes  or  if  the  parameters  change  drasticedly. 


3.2  Eigenvalue  and  Eigenfunction 

This  task  is  identical  with  itask=l  except  the  coefficients  of  the  Chebyshev  series 
and  the  values  of  the  variables  at  the  collocation  points  axe  evaluated.  Except  for 
itask,  the  specifications  axe  the  same  as  in  the  previous  section: 

(a)  One  line  with  itask: 

2  to  request  an  eigenvalue 

(b)  One  line  with  the  following  data: 

ievl  the  index  of  the  first  eigenvalue  parameter 

iev2  the  index  of  the  second  eigenvcdue  parameter  (for  complex 

problems  only) 

-1  an  optional  specification  to  reset  the  eigenfunction 

(c)  npar  lines,  where  each  line  specifies: 

par(k)  the  value  of  parameter  Pk-  The  values  of  par(ievl)  and 
par(iev2)  are  set  to  the  initial  guess. 
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The  real  and  imaginary  parts  of  eigenfunction  are  stored  in  file  0. function.  This 
file  may  be  use  by  PSH  as  initial  data.  The  eigenvector  is  normaJized  to  a  reasonable 
sccile  but  determined  only  within  a  complex  factor.  Amplitude  and  phase  of  eigen 
function  are  stored  in  O.famphafor  ploting  purposes. 


3.3  Tables  of  Eigenvalues 

This  task  enables  calculation  of  one-,  two-,  or  three-dimensional  tables  of  eigenvad- 
ues.  The  local  procedure  (it ask®  1,2)  must  be  performed  at  least  once  to  provide  a 
starting  point  for  this  task.  The  values  of  ievl  and  iev2  that  specify  the  eigenvalue 
are  taken  from  the  previous  task.  The  required  input  is: 


(a) 

One  line  with  it  ask: 

3 

to  request  a  table  of  eigenvalues  coefficients 

(b) 

One  line  specifying  the  dimension  of  the  table: 

ndim 

the  dimension  of  the  table  (integer) 

(c) 

ndim  lines. 

where  each  line  specifies: 

kpar 

the  index  of  the  parameter  to  be  varied  (integer) 

nstep 

the  number  of  steps  to  be  taken  (integer) 

step 

the  step  size  (real) 

The  first  parameter  specified  is  varied  in  the  innermost  loop.  The  eigenvalues  calcu¬ 
lated  are  stored  in  file  0. table. 


3.4  Curves  in  Parameter  Planes 

This  task  enables  the  direct  calculation  of  neutral  curves,  curves  of  constant  amplifi¬ 
cation,  and  other  useful  curves  in  the  plane  of  two  parameters  with  a  third  parameter 
completing  the  complex  eigenvalue  problem.  The  local  procedure  (itask=l,2)  must 
be  performed  at  least  once  to  provide  a  starting  point  for  this  task. 

To  describe  the  input  for  a  complex  problem,  we  imagine  a  three-dimensional 
coordinate  system  with  the  horizontal,  vertical,  and  normal  axes  formed  by  the  three 
parameters  pax(ivl) ,  par  (iv2) ,  and  par  (iv3) ,  respectively,  and  a  curve  in  the  plane 
spanned  by  par(ivl)  and  par(iv2).  Beginning  at  some  starting  point  provided  by 
the  local  procedure,  we  want  to  proceed  in  one  of  the  two  possible  directions  and  find 
new  points  of  the  curve.  All  but  the  three  parameters  used  as  coordinates  maintain 
the  starting  values  along  the  curve.  In  particular,  for  a  neutral  curve,  the  growth  rate 
in  the  starting  point  must  be  zero. 

Since  the  parameters  involved  may  have  different  orders  of  magnitude,  we  have  to 
choose  proper  scale  values.  For  this  step,  we  envision  a  plot  of  the  curve  on  a  sheet 
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of  paper  and  choose  the  scales  according  to  the  rules  of  nomography  so  that  the  finsJ 
curve  would  look  “nice”  without  being  squeezed  in  either  axis.  In  the  scaled  plot,  the 
distance  between  points  should  be  neither  too  large  nor  unnecessarily  small.  Also,  the 
plot  has  a  “window”  defined  by  the  minimum  and  maximum  of  the  variables  plotted. 
The  input  required  is  as  follows: 

(a)  One  line  with  itask: 

4  to  request  a  curve  in  a  parameter  plane 

(b)  Two  lines  specifying  data  for  the  “horizontal”  auid  “vertical”  direction,  respec¬ 
tively: 

ivl  ,2  the  index  of  the  associated  parameter  (integer) 
vgridl  ,2  the  scale  value  for  this  parameter  (real) 
vminl ,  2  the  minimum  value  (real) 
vmaxl  ,2  the  maximum  value  (real) 

(c)  For  complex  problems  only:  one  line  specifying 

iv3  the  index  of  the  second  eigenvalue  parameter  (integer) 

(d)  One  line  specifying  direction  and  step  size  along  the  curve  in  scaled  variables: 
angle  the  initial  direction  in  the  plane  of  the  curve  as  the  positive 

(counter-clockwise)  or  negative  angle  measured  in  degrees 
from  the  horizontal  axis 

radius  the  distance  between  points  along  the  curve 

(e)  One  line  specifying  the  number  of  steps: 

npoints  the  number  of  points,  or  arc  length  in  multiples  of  the 
radius,  to  be  traced  along  the  curve 

The  eigenvalues  calculated  are  stored  in  file  0. curve. 

The  iteration  for  the  eigensolution  can  neither  use  par  (ivl)  (if  the  curve  tangent 
is  horizontal)  nor  par  (iv2)  (if  the  tangent  is  vertical)  as  the  first  eigenvalue  parameter 
but  internally  iterates  normal  to  the  curve  in  the  plane  of  par  (ivl),  par(iv2).  for 
complex  problems,  the  second  eigenvalue  parameter  is  par(iv3). 

A  safe  recipe  for  the  choice  of  the  various  data  cannot  be  given  here  since  it  depends 
on  the  unknown  properties  of  the  curve.  A  reasonable  choice  for  the  scale  values 
would  be  the  order  of  magnitude  or  the  physical  values  of  the  parameters  par  (ivl), 
paT(iv2)  of  some  point  at  the  curve.  For  plane  Poiseuille  flow,  with  pair ( ivl )=Re, 
par(iv2)=  Or,  and  par(iv3)=  Uri  the  choice  vgridl=  10000,  vgrid2=  1  would  be 
almost  the  same  as  vgrid=  20000,  vgrid2=  2,  and  be  as  good  as  vgridl=  5772, 
vgrid2=  1.  Note,  however,  that  the  radius  is  the  distance  in  scaled  variables.  With 
a  radius  of  0.01,  the  maximum  physical  steps  in  the  horizontal  direction  would  be 
100,  200,  or  57.72  in  the  three  examples,  while  the  maximum  vertical  steps  would  be 
0.01,  0.02,  or  0.01.  The  steps  should  be  small  enough  to  obtain  a  smooth  plot  (and 
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rapid  convergence)  and  because  problems  may  arise  if  the  distance  between  multiple 
branches  of  the  curve  is  smaller  than  the  radius.  Neutral  curves  can  have  cusp-like 
sharp  turns,  as  for  the  boundary- layer  over  a  rotating  disk.  If  the  steps  2ire  too 
large  for  convergence,  or  if  the  specified  radius  exceeds  7r/8th  of  the  curve  radius  of 
curvature,  the  radius  will  be  automatically  divided  by  increasing  powers  of  2,  and 
hence  the  code  may  produce  2,  4,  or  more  points  per  step.  The  procedure  attempts 
to  double  the  reduced  internal  step  size  after  every  point  until  it  proceeds  with  the 
originad  radius.  Hence,  the  continuation  procedure  will  slow  down  in  regions  of  strong 
curvature.  This  internal  testing  consumes  additional  CPU  time,  and  one  may  as  well 
run  the  whole  curve  with  a  smaller  radius  instead. 

Ideally,  the  angle  (in  degrees)  would  be  the  positive  or  negative  angle  between  the 
horizontal  and  the  curve  tangent  in  the  starting  point.  The  two  signs  correspond  to 
the  two  directions  in  which  the  procedure  can  proceed.  The  choice  of  the  angle  is  not 
critical  as  long  as  the  radius  is  sufficiently  small  and  the  initial  direction  is  not  normal 
to  the  curve.  Trial  and  error  with  npoints=4  and  a  small  radius  will  help  (the  angle 
and  radius  for  the  last  point  are  printed).  Note  that  the  angle  changes  with  the  ratio 
vgridl/vgrid2. 

The  procedure  terminates  if  the  number  of  (full)  steps  exceeds  npoints,  if  the 
distance  between  any  point  (after  the  third)  and  starting  point  is  less  than  the  radius 
(to  prevent  unnecessary  loops  through  closed  curves),  and  if  the  curve  crosses  the 
window  given  by  the  minimum  and  maximum  values.  It  is  often  desirable  to  specify 
extreme  values  for  npoints  or  the  window  size. 


3,5  Close  and  Reopen  the  Input  File 

This  task  enables  to  suspend  the  program  execution,  edit  or  replace  the  Parameters, 
and  resume  execution  either  on  a  multiwindow  workstation  or  on  a  terminal  unaer 
Unix.  The  input  is  as  follows: 

(a)  One  line  with  it  ask: 

-1  to  request  a  pause  for  editing  the  Parameter  file 

The  program  will  close  the  file  Parameters  and  pause  until  continuation  is  requested. 
During  this  time,  the  file  can  be  changed.  After  continuation,  the  new  Parameters 
will  be  read  starting  at  the  top.  This  task  is  especially  helpful  in  the  initial  study  of 
a  new  proLlen  without  guidance  on  proper  parameter  ranges.  The  Parameters  file 
can  be  displayed  and  edited  in  one  window  while  test  runs  are  performed  in  another. 
Note  that  all  input  is  saved  in  the  input-log  file  0. input  for  later  inspection  and 
reruns. 
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3.6  Execute  System  Command 

This  task  enables  to  execute  system  commands  during  execution  of  the  code.  The 
input  is  as  follows: 

(a)  One  line  with  it  ask: 

-2  to  request  execution  of  a  system  command 

text  the  text  passed  to  the  system  Ctill 

Deviating  from  the  general  format,  the  text  may  contain  spaces  and  is  terminated  by  a 
slash  /  or  the  end  of  the  line.  The  commzind  is  typically  used  to  copy  or  move  output 
files  into  safe  files  before  they  are  overwritten  by  one  of  the  consecutive  tasks.  A 
typical  application  is  securing  part  of  a  curve  before  the  continuation  in  the  opposite 
direction  is  requested. 
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Chapter  4 

Control  for  PSH 

The  file  Control  is  read  in  a  format  similar  to  the  free  format  in  Fortran  77.  Excep¬ 
tions  are:  all  items  must  be  separated  by  spaces  (no  commas  are  allowed);  the  items 
may  be  character  strings;  and  the  data  cannot  be  continued  on  the  next  line.  The 
slash  /  terminates  the  reading  of  data  from  the  line.  Empty  lines  axe  ignored  and  can 
be  freely  used  to  enhance  readability. 

The  file  Control  consist  of  various  groups  of  data  described  in  the  following: 

(1)  A  single  line  of  text  specifying  the  problem. 

(2)  One  line  specifying  the  following  four  data 

nstart  station  number  at  which  initial  data  is  specified 

nf  inal  final  station  number 

ndelta  station  numbers  per  step 

step  integration  step  size  in 

(3)  One  line  specifying  the  integer 

inlt  the  flag  controlling  whether  the  nonlinear  terms  are  implicitly 
evaluated  at  the  current  step  or  explicitly  at  the  previous  step 

0  explicit 

1  implicit 

(4)  One  line  specifying  the  integer 

itx  the  maximum  number  of  iterations  to  be  performed  for  each 
mode  per  step 

(4)  One  line  specifying  the  following  two  data 

tpar  temporal  frequency  of  the  fundamental  mode 
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zpar  wave  number  in  the  coordinate  direction  for  the  fundamental 

mode 

One  line  specifying  the  two  indices 

itr  temporal  index  of  the  reference  mode 

izr  index  of  the  wave  number  in  the  (f^  direction  for  the  reference 

mode 

One  line  specifying  the  following  two  thresholds 
thron  threshold  maximum  magnitude  of  the  nonlinear  forcing  terms 
for  a  new  mode,  above  which  a  new  mode  is  created  and 
integrated 

throf f  threshold  amplitude  of  the  normv  disturbance  variable  below 
which  an  existing  decaying  mode  will  be  discarded 

One  line  specifying  the  number  of  mode-control  lines  nmgr. 

nmgr  lines,  where  each  line  specifies  control  of  a  mode  by  the  following  data: 
it  the  temporal  frequency  index  of  the  mode  in  multiples  of  tpar. 

iz  the  ^^-wavenumber  index  of  the  mode  in  multiples  of  zpar. 

ngr  mode  control  parameter  where 

0  indicates  that  the  mode  is  calculated  by  including  non¬ 
linear  effects  of  other  modes  on  itself  but  neglecting  its 
effect  on  other  modes 

1  the  mode  is  calculated  and  its  interaction  with  other 
modes  is  included,  and  the  real  part  of  its  wavenumbers 
is  determined  by  the  pheise-locking  condition  with  the 
reference  mode 

2  the  mode  is  calculated  and  its  interaction  with  other 
modes  is  included,  and  the  real  part  of  its  wavenumber 
is  determined  independently  by  imposing  a  restriction 
on  its  norm  (default) 

3  the  mode  is  calculated  and  the  nonlinear  effects  on  this 
mode  is  neglected  by  setting  the  right  hand  side  to  zero. 

Its  wavenumber  is  determined  by  imposing  a  restriction 
on  its  norm 

-1  indicates  that  the  mode  is  to  be  neglected  altogether 

One  line  specifying  the  number  of  initial  modes  nmodes. 

nmodes  lines,  where  each  line  specifies  an  initial  mode  by  the  following  data: 
it  the  temporal  frequency  index  of  the  mode 


iz  the  ^^-wavenumber  index  of  the  mode 

Digs  mode  control  option  for  the  initial  mode,  as  described  above 
in  item  7.  This  value  of  mgs  will  overwrite  mgr  for  this  mode 
if  it  is  specified  in  item  7. 

ampO  amplitude  of  the  initial  mode  for  the  mvar  variable  chosen  for 
imposition  of  the  norm  constraint 

phsO  phase  of  the  initial  mode  for  the  normv  variable  chosen  for 
imposition  of  the  norm  constraint 

Not  that  the  indices  of  the  modes  specified  here  should  be  within  the  limits  specified 
by  parameters  (modt)  and  (modzm,  modzp)  specified  in  common. ins.  A  summary  of 
output  data  for  every  step  of  psh  rtm  is  written  out  to  the  file  0 .  summary. 
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Chapter  5 
Installation 


The  code  is  delivered  in  tar  format  on  a  cartridge  tape  compatible  with  Silicon  Graph¬ 
ics  standard  cartridge  tape  drives.  To  provide  full  functionality,  the  code  should  be 
installed  in  the  user  directory.  No  root  privileges  axe  required.  Insert  the  tape  car¬ 
tridge  into  the  drive.  To  conform  with  the  original,  make  a  new  directory  w.psh  amd 
read  the  tar  file: 

$  mkdir  w.psh 
$  cd  w.psh 

$  tar  xvo  (  or  tax  xvof  /dev/???  ) 

$  Is 

If  the  default  device  /dev/tape  is  not  available  on  your  workstation,  use  the  proper 
device  name  for  the  cartridge  tape  drive  on  a  remote  workstation.  After  the  instal¬ 
lation  is  complete,  there  will  be  one  file  README,  a  Makefile,  a  compressed  tar  file  of 
the  code  that  should  be  saved,  amd  a  series  of  directories.  First,  read  README  to  get 
an  update  on  chamges,  bugs,  or  other  news. 

To  test  for  successful  installation,  execute  the  following  commamds: 

$  cd  c.cone 
$  Is 
$  maike 
$  Ish 
$  Is 

The  compairison  of  the  file  lists  should  show  a  series  of  output  files  such  as  0 .  input. 
The  names  of  all  output  files  staurt  with  0 .  which  simplifies  their  handling  or  removal. 
If  the  code  ran  error-free  to  completion,  the  output  files  should  be  identical  with  those 
in  the  subdirectory  w. results,  except  for  the  date.  To  verify  proper  function  of  the 
code,  check  with 
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$  diff  .  V. results 

This  command  should  reveal  no  diiferences  for  0.*  files.  The  command 
$  make  clean 

removes  the  problem  specific  object  files  in  ../w.objs,  the  executable  files  Ish  and 
psh,  and  prompts  whether  you  want  to  delete  the  output  files  in  the  current  directory. 
New  applications  must  be  prepared  in  directories  on  the  same  level,  i.e.,  in  directories 
w.psh/c.new-appl  to  be  compatible  with  the  various  Makefiles. 
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Chapter  6 

Examples  of  Applications 


6.1  Hypersonic  Flow  over  a  Blunt  Cone 

Insert  files  and  input  files  for  studying  stability  of  hypersonic  flow  over  a  blunt  cone 
are  set  up  in  directory  c.bluntcl.  These  files  are  used  as  examples  in  this  section  to 
illustrate  how  the  user  may  build  his  own  physics  insert  files  and  input  files  to  tackle 
a  new  problem. 

6.1.1  metric.ins 

Metric  terms  for  transformation  from  Cartesian  coordinates  (x,y,  z)  to  axisymmetric 
body-intrinsic  coordinates  =  s,  =  rj,  =  <!>)  and  local  Cartesian  coordinates 
aligned  with  the  cone  surface  are  given  here.  They  are  coded  into  the  metric.ins 
insert  file. 


Body-Intrinsic  Coordinates 

This  coordinate  system  (Figure  6.1)  is  used  when  icoord*!  is  specified  in  Definitions. 
Transformation  is  from  Cartesian  coordinates: 


to: 


=  5;  arc- length  along  the  body 
=  Tj;  normal  distance  from  the  body  surface 
=  <f>;  azimuthal  angle 
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Figure  6.1:  Body- Intrinsic  Coordinate  System  Used  for  Stability  Analysis 


X  =  --  risin{0c{s))} 

y  =  +  7Cos(^c(5))}  sin{<i>) 

z  =  {'’6(5)  +  »?COs(^c(5))}  COs{(f>) 

where  rj,(s)  is  the  radius  of  the  body,  Xj,(s)  is  the  axial  distance,  and  tan{0c('S)}  is 
the  slope,  at  the  body  surface  (t/  =  0). 

First  derivatives  are  : 

dx 


dxb 

=  —-r)cos{0As\)  — 


dx^ 

w 


dx^  _  dy 

~  Ts 

=  -  V  sin{Oc{s))  sin{<p) 

dx^  _  dy 
drj 

=  cos(^c(s))5zn((^) 

dx"^  _  dy 

le  ~  ^ 

=  {rfc  +  j7co5(^c(5))}  cos{4>) 


dx^ 

W 


dx^ 

W 

dx^ 

W 

Second  derivatives  axe; 


dz 

d^ 

-  »/ stn(^c(s))  ^}cos{4>) 
dy 

cos(^c(5))  cos{(f) 

dz 

d4> 

-  {n  +  7?cos(^c(s))}  $in{<^) 


d^x^  _  d^x 
d(^  ~  ^ 

=  +  n  sin{ec{s)) 

d^x^  _  d^x 
di'^d^^  dsdy 


,  *  /  \  \  r 

=  -  COs(^c(5)) 


d^x^ 

'  '  "  ds 
d^x 

d^de 

dsd(l> 

=  0 

d'^x^ 

d'^y 

di^dV- 

ds^ 

S^Ii 


ycos{9^{s)) 


—  T]COs{9c{s)) 


<P9c 

ds^ 


"  1^1  ■  * 
-  T}  stn{9c{s))  ^ 
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d^de 


d^x^ 

d^de 


d^y 

dsBrj 


—  sin{Oc{s))  sin  4> 
as 

d^y 


dsd4> 

3in(de(5))  —}  cos{(i>) 

d'^y 

dr}d<i> 

cos{dc{s))  cos  4> 
d'^y 
d<t>d(i> 

-{n  +  T]  cos{6c{s)) }  sin{<f>) 


d^x^ 

d^de 

d^x^ 

dede 

d^x^ 

d^x^ 

5^35^3 


ds^ 

{^-Vcos{ec{s)) 

d^z 

dsdrj 

—  sin 

dh 


(de.V 

(irj 


—  sin{Oe{s))  cos  <l> 
as 


dsd<i) 

-{^  -  ?/ sm(5c(5))  sin{<f>) 

d'^z 

dr}d<f> 

—  cos{Oc{s))  sin  <f> 
d^z 
d<t>d<t> 

-{vb  +  Ti  cos(6c{s)) }  cos{<i>) 


For  the  straight  cone  region,  with  uniform  grid  distribution  in  s,  there  axe  only 
two  nonzero  third  derivative  terms. 

d^x^  d^y 

5^15^35^3  -  dsd<i>d<t> 

=  -vsin{ec{s))^}sin{<i>) 
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5^X2  diy 

~  drfd4>d<l> 

=  —  cos{6c{s))  sin  <f> 


Local  Cartesian  Coordinates 

This  coordinate  system  is  used  for  stability  analysis  when  icoord=0  is  specified.  Here, 
all  curvature  terms  and  cone  divergence  are  neglected. 

X  =  s 

y  =v 

z  =  z 

dx^  dx  _ 

—  _  ^  , 

de  ~  dri~ 

dx^  _  _  1 

W  ~ 

Since  the  basic  flow  is  axis5mametric,  and  disturbamces  are  periodic  in  direction 
with  wave  number  n^,  ^  is  set  to  zero,  with  <^  =  0  in  the  above  metric  terms. 

6.1.2  gridbs.ins 

Subroutine  gridbs  assigns  the  basic  state  grid  simply  by  reading  it  in  along  with 
the  basic  state  solution.  Auxiliary  routines  that  gridbs  calls  to  accomplish  this  are 
described  more  fully  in  Appendix  C.  Of  these,  three  axe  specific  to  the  basic  state 
flow  solver  used.  They  are;  bsiois,  bpstda,  and  scale  . 

Subroutines  bsiois  amd  bpstda  read  in  the  station  data  and  flow  profiles,  respec¬ 
tively,  from  a  FORTRAN  formatted  file  bpse .  inp.  This  data  file  contains  the  basic 
state  solution  amd  derivatives  generated  by  the  Thin  Layer  Navier-Stokes  (TNLS) 
code  (Esfahanian  1991).  The  TNLS  code  uses  a  nonorihogonal  axisymmetric  body- 
intrinsic  grid.  The  equations  axe  solved  in  the  computational  domain  (I,J).  The 
outer-most  grid  line  in  the  streamwise  direction  (^{JMAX)  is  the  bow  shock.  Grid 
lines  in  the  streaunwise  direction  (J-|-constaint)  aure  generated  by  dividing  the  distance 
between  the  shock  and  the  cone  surface  by  a  stretching  function  (stored  in  am  array 
an(J)).  Grid  lines  (I=constant)  axe  normal  to  the  cone  surface. 
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bsiois 


Subroutine  bsiois  reads  in  station  data  and  profile  data  and  stores  the  station  data  in 
the  array  stndat.  Profile  data  are  read  from  the  sequential  data  file  for  error  checking. 
It  also  sets  the  scales  for  the  basic  state  solution  by  calling  the  subroutine  therm  (see 
appendix).  The  relevant  open,  read  and  station  data  assignment  statements  are 
given  below. 

i  *  inpprf 

open(unit=i,f ile=profil,access®*sequential’ .form® ’formatted’ , 

I  status® ’old’) 
revind(i) 

c  read  in  all  the  station  data  (except  the  profiles  themse.  /es)  from 
c  input  and  store  them  . . .  (assime  that  station  number  nstn-*-! 
c  follows  station  number  nstn  in  the  input  file.) 

lcp=len(charv) 
do  10  iv®l,nc 
do  10  k®l,lc 

write(charv(iv) (k:k) , ’ (al) *)  *  ’ 

10  continue 

read(i,’(a80)’)  charv(l) 
c  tcone  ®  cone  half  angle  (degrees) 

c  scone  ®  arc  length  along  sxirface  at  end  of  spherical  nose  (m) 
c  aminf  ®  free-stream  mach  niunber 
c  reinf  ®  free-stream  re3rnolds  number 
c  tinf  ®  free-stream  static  temperature  (k) 

c  pinf  ®  free-stream  static  pressure  (psia) 

c  reflen  ®  spherical  nose  radius  (m)  (reference  length) 
c  betaa  ®  basic  state  grid  stretching  paurameter 

c  tw  ®  basic  state  grid  stretching  parameter 

c  pr  ®  prandtl  number 

read(i,*)  tcone .scone .aminf .reinf .tinf .pinf .reflen. 

I  betaa. tw.pr 

c  specify  the  free-stream  mach  number,  reference  length,  and  inverse 
c  reynolds  number  by  making  the  I’s  digit  in  mdlair  equal  to  1  and  the 
c  10 ’ s  digit  equed  to  1  ... 

mdlair  =  iabs(nint(con(2))) 

mdlair  ®  mdlair  -  mod (mdlair. 100)  +11 


c  compute  the  ratio  of  specific  heats  (gamma),  reference  velocity,  ... 
c  set  the  local  pressure  equal  to  1  and  disregard  the  equation  of  state 
c  data  here . 

fsrei  -  one/reinf 

call  therm  (mdlair,0,tinf ,fspres,bsrrho,bsrvel,reflen,dypdep, 

I  aminf , f srei ,pr ,gae , cpe .dyvisc ,sve , coe .one .one ,rOs , 

I  r Is , r2s , cpO , cp 1 , cp2 , cp3 , cp4 , f vO . f V 1 , f v2 , f v3 , f v4 , 

I  svO , svl , sv2 , sv3 , sv4 , coO , col , co2 , co3 , co4) 


c  vahid’s  velocity  scale  is  sqrt(r*tinf ) ,  where  r  *  gas  constant. 

bsrvel  =  bsrvel/sqrtCgae) /aminf 

realv(  1)  *  reflen 
reailvC  2)  =  reflen 

c  xi“-C3}  =  circumferential  angle  (dimensionless)  . . . 

realv(  3)  =  one 

re2dv(  4)  =  reflen/bsrvel 

realv(  5)  =  bsrrho 
realv(  6)  =  tinf 
realv(  7)  =  bsrvel 
realv(  8)  =  bsrvel 

realv(  9)  *  pinf 
realv(lO)  -  gae 
realv(ll)  =  aminf 
realv(12)  =  reinf 
re€dv(13)  *  pr 
realv(14)  =  tcone 
readvdS)  -  scone 
realv(16)  =  betaa 
realv(17)  =  tw 

j  =  0 

17  j  =  j  +1 

if (j+l.le. jmax)  j=j+l 

read(i,902,iostat=ie(l))  an(j) ,dadc2(j) ,d2adc2(j) ,d3adc2(j) 
go  to  17 

902  format(lx,4e23.16) 
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c  store  the  stretching  function  and  its  derivatives  in  realv  . . . 


nr0»18 

do  20  j*0,ja 
realv (nrO+j) 
realv (nr0+ j  + j  a+ 1 ) 
realv (nr0+ j  +2* j  a+ 1) 
readv (nr0+ j  +3* j  a+ 1 ) 
20  continue 


an(j) 

dadc2(j) 

d2adc2(j) 

d3adc2(j) 


nrp*17+4*(ja+l) 


c  read  to  the  next  station  while  checking  for  errors  in  the  input  file 
istn  *  0 
nsloc  >  1 
SO  continue 

read(i,901 ,iostat=ie(l) , end-900) 

I  nstn,np,arclen, 

c  dependent  variables  at  shock 

I  qs,rhos,ts, 

c  grid  data  at  body 

I  xb , yb , xbdc 1 , ybdc 1 , xbdc2 , ybdc2 , xbdc3 , ybdc3 , 

c  grid  data  at  shock 

I  xs , ys , xsdc 1 , ysdc 1 , xsdc2 , ysdc2 , xsdc3 , ysdc3 


do  55  n^O.np 

read(i,911,iostat*ie(l))  nread.xij ,yij , 

1  bsoln(l) ,dbdc(2,l) ,d2bdc2(2,2,l) , 

I  bsoln(2) ,dbdc(2,2) ,d2bdc2(2,2,2) , 

I  bsoln(3) ,dbdc(2,3) ,d2bdc2(2,2,3) , 

I  bsoln(4) ,dbdc(2.4) ,d2bdc2(2,2,4) 

55  continue 

c  if  the  entire  dataset  has  been  read  successfully,  increment  the  loop 
c  counter  and  store  the  station  data  . . . 

istn»istn+l 
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c  shock  standoff  distance 

sns  ®  dsqrt((xs-xb)**2  +  (ys  -yb)**2) 

stndatC  0,istn)  ^  nsloc 

stndatC  l,istn)  -  arclen 

StndatC  2,istn)  =  np+1 

StndatC  3,istn)  =  0  !  xi‘{3} 

StndatC  4,istn)  =  reflen 

StndatC  5,istn)  =  qs 

StndatC  6,istn)  >  rhos 

StndatC  7,istn)  *  ts 

StndatC  8,istn)  -  sns 

StndatC  9,istn)  =  xb 

StndatC 10, istn)  =  yb 

StndatC 11. istn)  =  xs 

StndatC 12, istn)  =  ys 

StndatC 13, istn)  =  xbdcl 

StndatC 14, istn)  *  ybdcl 

StndatC IS, istn)  =  xbdc2 

StndatC 16, istn)  *  ybdc2 

StndatC 17, istn)  =  xbdc3 

StndatC 18, istn)  *  ybdc3 

StndatC 19 « istn)  -  xsdcl 
stndatC20,istn)  =  ysdcl 
stndatC21,istn)  =  xsdc2 
StndatC 22, istn)  =  ysdc2 
stndatC23,istn)  =  xsdc3 
stndatC24,istn)  =  ysdc3 

c  end  if 

nsloc  =  nsloc  +  Cnp+1) 
go  to  50 

900  continue 

c  close  the  input  file  . . . 
close Ci) 


901  formatC2i5,20e23.15) 
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911 


f oxmat (ilO , 14e23 . IS) 


bpstda 

This  subroutine  opens  the  sequential  data  set  named  profil,  connects  it  to  unit 
inpprf ,  and  reads  the  basic  state  profiles.  It  tramsforms  the  Cartesian  velocity  com¬ 
ponents  calculated  by  the  TNLS  code  to  physical  contravariant  components  in  the 
stability  analysis  coordinates.  Then  it  writes  out  these  profiles  to  the  direct  access 
file  named  daprof  which  has  already  been  opened  and  connected  to  imit  iscr  by 
bdiois.  The  relevant  open  and  read/write  statements  are  given  below. 


logical* (llogic)  logicvCO:*) 
chaa-acter  chaiviO :*)*{*) 
integer* (lint)  intv(-6:*) 

real*(lreaLl)  stndat(0:nstapr.inznsta) ,realv(0:*) 
dimension  ie(*) 

c  nts  -  ntimber  of  theTmod3rnafflic  state  variables  in  equation  of  state 
c  nbsv  =  ntnaber  of  basic  state  variables, 
c  j2  «  number  of  basic  state  grid  points. 

parameter  (one»l,2ros0,nts=2,nbsv=4,j 2=199) 

reaQ.*Clreal)  an,z,bsoln(nbsv) ,dbdc(4,nbsv) ,d2bdc2(4,4,nbsv) 

chauracter  profil*(*) ,daprof*(*) 

i  =  inpprf 

open (unit=i,file=profil, access* 'sequent iad’ .form* 'formatted' , 

I  status* 'old') 
revindCi) 

read(i, ’ (aSO) ’)  charv(l) 

read ( i , * )  t cone , scone , aminf , reinf , tinf .pinf , ref len . 

I  betaa.tv.pr 

c  wall  normal  stretching  function  and  its  derivatives  .... 

j  =0 

ja=0 

read(i.902,iostat*ie(l))  an 
902  format(lz,4e23.16) 

anl  gATi 
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17  if  (ie(l).eq.O)  then 

write (6, 902)  an 

if  (abs(one-an) .le.l.e-3*abs(an-anl))  go  to  18 
else 

write(6,  ’  (//,lh  ,72(lh*)  ,//,3x. ’’subroutine  bpstda:  error”, 

I  ”  occurred  while  reading  point  j*”  ,i3,/,3x,  ”of  ”, 

I  ’’stretching  function.  stop!”,//,lh  ,72(lh*) ,//)  ’ )  j 

stop 
end  if 
anl^an 
ja=ja+l 

if  (j+l.le.j2)  j=j+l 
read(i,902,iostat=ie(l))  an 
go  to  17 

18  if  (j.ne.ja)  then 

write(6,  ’  (//,lh  ,72(lh*)  ,//,3x, ’’subroutine  bpstda:  increase”, 
I  ’’  dimension  j2  to  at  least  ’’,i3,’’.  stop!”,//, 

I  Ih  .72(lh*))’)  ja 

stop 
end  if 

c  number  of  stations  in  file  . . . 

read(i,909,iostatsie(l))  nstmax 
909  format (ilO) 

if  (ie(l).ne.O)  then 

write(6, ’ (//,lh  ,72(lh*) ,//,3x, ’’subroutine  bpstda:  error  ”, 

I  ’’occurred  while  reading  total  number”  ,/,3x,  ”  of  ”, 

I  ’’stations.  stop!”,//,lh  ,72(lh*) ,//) ’) 

stop 
end  if 

c  read  in  all  the  station  data  (except  the  profiles  themselves)  from 
c  input  and  store  them  . . .  (assume  that  station  number  nstn+1 
c  follows  station  number  nstn  in  the  input  file.) 

pi  -  4*datan(l.d0) 
thetac  *  tcone*pi/180.d0 

gaonma  -  reaQ.v(10) 
aminf  =  realv(ll) 

do  300  k*l,intv(2) 
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readCi, ’ (2iS) ’)  nstn.np 
key  s  nint(8tndat(0,k)) 
do  8S  n«0,np 

do  80  iS^l.nbsv 
bsoln(i3)  «  zro 
do  75  i2*l,4 
dbdc(i2,i3)  »  zro 
do  70  il*1.4 
d2bdc2(il,i2,i3)  »  zro 
70  continue 
75  continue 
80  continue 

read(i,911,iostatsie(l))  nread,zij ,yij , 

bsoln(l) ,dbdc(2,l) ,d2bdc2(2,2,l) , 
bsoln(2) ,dbdc(2,2) ,d2bdc2(2,2.2) , 
bsoln(3) .dbdc(2.3) ,d2bdc2(2,2.3) , 
bsoln(4) ,dbdc(2.4) ,d2bdc2(2.2.4) 

format(il0,14e23. 15) 

if  (ie(l).eq.O)  then 

c  assume  for  nou  that  the  angle  thetac  between  the  body  axis  and  the 
c  surface  is  a  constant. 

if  (n.eq.O)  rb*yij 
z  =  (yij-rb) /cos (thetac) 

c  transform  the  velocity  components  from  aziad  and  radial  to 
c  tangential  and  normed. 

utan  =  bsoln (3)  *dcos  (thetac)  bsoln(4)*dsin(thetac) 
unor  =  -bsoln(3)*dsin(thetac)  +  bsoln(4)*dcos (thetac) 

bsoln(3)  =  utan 
bsoln(4)  =  unor 

write(iscr ,rec*key+n)  z, (bsoln (i3) ,i3=l ,nbsv) 
else 

c  shoTildn’t  be  any  errors  because  subroutine  bsiois  already 
c  successfully  read  in  file. 

write(6,*)  ’stop:  error  in  bpstda.’ 


I 

I 

I 

I 

911 
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stop 
end  if 

nstnl  -  nstn 

85  continue 

300  continue 

return 

end 


6.1.3  bstate.ins 

Subroutine  bstate  reads  in  and  interpolates  the  basic  state  profiles  from  the  direct 
access  file  created  by  subroutine  gridbs.  The  method  by  which  it  does  this  is  inde¬ 
pendent  of  the  input  and  output  routines  required  to  read  in  the  basic  state,  and  so 
need  not  be  changed  for  different  basic  state  flow  solvers.  It  is  described  more  fully 
in  Appendix  C. 

6.1.4  btrans.ins 

This  insert  file  is  written  to  transform  basic  state  velocity  components  specified  in 
bstate.ins  into  either  Cartesian  or  physical  contravariant  velocity  components  (of  the 
stability  coordinate),  and  their  derivatives  into  derivatives  with  respect  to  stability 
coordinates.  As  physical  contravariant  velocity  components  of  the  basic  flow  are 
defined  here,  IBVEL  =  1  is  specified.  If  the  basic  flow  is  assumed  to  be  parallel 
(mcoord»0)  then  normal  (to  the  surface)  component  of  velocity  and  its  derivatives 
with  respect  to  t]  are  set  to  zero.  All  the  streamwise  derivatives  (with  respect  to  s) 
are  also  set  to  zero.  If  mcoord=l  the  aforementioned  nonparallel  terms  of  the  basic 
state  are  included  in  the  analysis. 

Also  defined  in  this  insert  file  are  the  equation  of  state  for  a  perfect  gas 


and  specific  enthalpy 


p  =  ^ 

■ 

h  =  CpT 


and  their  derivatives  with  respect  to  thermodynamic  variables.  Specific  heat  at  con¬ 
stant  pressure  Cp,  viscosities  /ui  and  H2i  and  thermal  conductivity  k  of  the  basic  state, 
and  their  derivatives  with  respect  to  temperature  are  specified  by  calling  therme  en¬ 
try  statement  in  subroutine  therm  (see  Appendix  B).  Constant  Prandtl  number  of 
0.72,  specific  heat  ratio  7  =  1.4,  and  the  Sutherland  Law  for  viscosities  are  used. 
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6.1.5  boundary.ins 

The  present  insert  file  is  set  up  so  that  the  following  sets  of  conditions  may  be  chosen 
by  the  user  by  specifying  the  axray  idbc(nvar)  in  the  input  file  Definitions. 

1.  homogeneous  conditions  at  both  boundaries  for  both  stationzury  and  oscillatory 
modes;  idbc  =  1. 

2.  homogeneous  conditions  at  both  boundairies  for  both  stationary  and  oscillatory 
modes  except  temperature  gradient  of  stationary  modes  are  set  to  zero  at  the 
wall;  idbc  =  2. 

3.  homogeneous  conditions  at  the  wall  for  both  stationary  and  oscillatory  modes 
and  asymptotic  boundary  conditions  at  the  faxfield;  idbc  =  3. 

4.  homogeneous  conditions  for  both  stationary  and  oscillatory  modes  at  the  wall 
except  for  normal  temperature  gradient  and  asymptotic  boundary  conditions  at 
the  farfield;  idbc  =  4. 

If  homogeneous  conditions  are  chosen,  the  outer  boundary  should  be  located  at  a 
larger  distance  from  the  wall  than  that  used  with  asymptotic  conditions. 

6.1.6  pointd.ins 

Most  of  the  results  for  this  flow  are  obtained  with  ffmax  varying  as  y/s.  The  local 
variable  btrf  for  Tjmas  used  in  subroutine  f  dinit  is  vairied  by  the  following  statement 
in  pointd.ins. 

btrf  =  strfp(2)*sqrt(  par(l)*par(2)  ) 

6.1.7  postpr.ins 

The  present  module  is  set  up  to  calculate  the  following  new  variables: 

1 .  pressure  disturbance,  rp  (n ,  it ,  iz) 

2.  massflow  disturbance,  rmass(n,it  ,iz) 

3.  total  temperature,  t  ot  (n ,  it ,  iz) 

4.  j  —  th  covariant  component  of  wall  shear  stress  for  ste2uiy  {it  =  0,  iz)modes, 
wshcovCj ,iz) 

5.  wall  heat  transfer  for  steady  {it  =  0,  iz)  modes,  in  j  —  th  covariant  direction, 
dqwalKj  ,iz) 

6.  m20cimum  amplitude  for  massflow  fluctuation  for  the  reference  mode  {itr,  izr) 
and  corresponding  growth  rates  of  other  modes  at  the  same  location,  ampmax  (it ,  iz) 
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7.  growth  rate  profiles  (n  =  0,  ji)  of  ig-ih.  variable  at  the  location  where  the  mass- 
flow  disturbance  of  the  reference  mode  {ttr,  izr)  is  maximum,  grate  (a ,  it ,  iz ,  ig) . 
Here  ig  =  1,..,4  representing  growth  rates  for  p,  T,  and  mass  flow  m,  re¬ 
spectively. 


6.1.8  Definitions 

A  S2unple  input  file  Definitions  for  stability  analysis  is  given  in  Figure  6.2.  Since 


stab .  eqs 
200 
10 
13 

-1  -1  - 
0  4  4 

-1 

4,  0.000, 
rbo  , 
T 

u  , 

V  , 

w 

WallProp 

Air-Pr«® 

Tinf ty 

Pinfty 

Lre£ 

Prandtl 

BlAodel 

Hflux 

tqpprt 

rela^v 

NGl 

NG2 

NG3 


L  rho,  T,  phy.  contr.  u  in  body- intrinsic  coordinates 
/no.  o£  collocation  points  -1 
/nvunber  o£  constants 
/nxinber  o£  parameters 

1  -1  -1  /idsy(nvar)  -1  eta  no  syasnetry,  0  sysm,  1  antisyocn 

4  4  /idbc(nvar)  type  o£  boundary  conditions 

/idpts^l  Gauss-Lab.,  0  b-state  grid,  -1  use  idtr£ 
22.,  1.12,  0.  /idtr£,  atr£,  btrf,  str3,  str4 
1  /name,  z-synmetry  o£  variable  1 
1  /name,  z-symmetry  o£  variable  2 
1  /name,  z-sysnetry  o£  variable  3 
1  /name,  z-symmetry  of  variable  4 
-1  /name,  z-syosnetry  of  variable  5 


0.0000000000000 
1120103. 
54.30000000000 
4.082743600e-03 
0.6667S0000e-^00 
0.7200000000000  /  con(6) 
311.00000000000 
O.OOOOOOOOOE-04 
5.000 
0.000 

6.000000000000 
6.000000000000 


/  con(l),  not  used 
/  con(2),  model  4  for  air  properties 
/  con (3),  Freestream  tetqpexature ,  (Kelvin). 

/  con(4),  Freestream  pressure,  (atm). 

/  con(5).  Reference  length  (m)  . 

Edge  Prandtl  ntimber. 

/  con (7),  Analy.  model;  invar,mcoord,mdlddx 
/  con(8).  Heat  flux  (not  used) 

/  con(9),  print  freq.  in  steps  for  profile  data 
/  con(lO),  1  dp/dx  removed  for  oscil.  modes,  0  NOT 
/  points  used  to  interpolate  basic  state  in  xi^d) 

/  points  used  to  interpolate  basic  state  in  xi^(2} 

6.000000000000  /  points  used  to  interpolate  basic  state  in  xi^{3} 


xi^'d) 

,  2, 

0 

/par(l) , 

affects 

Reinv 

,  0, 

0 

/par(2) , 

appears 

alphar 

,  1, 

1 

/par{3) , 

appears 

alpha! 

,  1, 

1 

/par  (4) , 

appears 

beta 

,  1, 

0 

/par(5) , 

appears 

omegar 

,  0, 

1 

/par(6) , 

appears 

omega! 

,  0, 

2 

/par (7) , 

appears 

UsqdRT 

,  2, 

0 

/par (8) , 

appears 

x!'{3) 

,  2, 

0 

/par(9) , 

affects 

complex (imag) 


Figure  6.2:  Example  Input  File  Definitions  for  Stability  Analysis  of  Hypersonic 
Flow  Over  a  Blimt  Cone  at  M^o  =  8- 


idbc  is  set  to  4,  asymptotic  boundary  conditions  axe  used  at  farfield.  Here  a  slight 
logarithmic  stretching  (idtrf=4  with  strfp(3)  =1.12)  is  used.  The  wall  is  at 
^Ttitn  =  strfp(l)  =  0,  while  the  farfield  boundary  at  the  reference  station  is  lo¬ 
cated  at  nondimensional  f/mox  =  strfp(2)  =  22.  Note  that  a  larger  value  should 
be  used  for  strfp(2)  if  homogeneous  conditions  are  applied  at  farfield.  In  the  next 
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nveur— 5  lines,  the  names  of  dependent  variables  for  disturbances  auad  their  symmetry 
property  in  directions  are  defined. 

Constants  and  parameters  may  be  defined  by  the  user  to  suit  a  particular  appli¬ 
cation.  There  axe  ncon=13  constants  and  npar=9  parameters,  con(l)  is  not  used, 
con (2)  named  as  Air-prop  specifies  the  thermodynamic  model  for  the  air.  This 
is  specified  as  1120103.  (see  Appendix  A  for  details).  Constants  con(3)  and  con(4) 
give  reference  temperature  and  pressure  for  nondimensionalization.  Their  edge  values 
vary  with  streamwise  distance  for  a  blunt  cone,  unlike  those  for  a  sharp  cone.  There¬ 
fore,  fixed  basic-state  freestream  values  of  temperature  (tinf  ty  =  =  54.3  K)  and 
pressure  (pinfty=p^  =  0.004082743  atm.)  Me  used  as  reference  values.  The  next 
constant  con  (5)  is  the  location  of  the  reference  station  in  arc  length  (meters)  and 
is  named  as  Lref.  In  our  results,  the  station  at  s’/Rff  =  175,  s*  =  0.66675m.  was 
input  as  a  reference  station.  Since  the  10°  digit  in  the  Air-Prop  model  is  set  to  3,  the 
reference  length  scale  rl,  is  the  boundary-layer  length  scale  =  2.8511  x  lO""*  m 
at  the  reference  station,  s*/Rs  =  175.  Prandtl  number  is  specified  as  con(6). 

The  three  digits  of  the  real  constant  blmodel  =con(7)  specify  the  following: 

10^  digit  mvar,  the  choice  of  dependent  disturbance  variable.  Dif¬ 
ferent  choices  are  described  in  Section  2.2.6.  The  meaning 
of  this  flag  mvar  is  generic  for  any  application. 

10^  digit  mcoord,  the  choice  of  coordinate  system  for  stability  anaJ- 
ysi' 

0  locail  Cartesian  coordinates  tangent  to  the  cone 
surface 

1  curvilinear  coordinates 

10°  digit  mddlx,  the  flag  for  basic-state  nonparallel  terms 

0  parallel  basic  state,  velocity  component  and 
gradients  of  baisic  flow  axe  neglected 
1  nonpairaJlel  basic  state,  velocity  component  and 

gradients  of  baisic  flow  au-e  included 

In  the  example,  physicaJ  contravaxiant  velocities  (mvau:=3)  in  curvilinear  coordinates 
(mcoord=l)  for  a  nonparallel  basic  state  (mddlx=l)  are  chosen.  The  rest  of  the  con¬ 
stants  axe  self-explanatory  by  way  of  comments  in  the  example. 

The  nine  dimensionless  parameters  defined  in  this  example  represent  the  following: 

par(l) 

pair  (2)  Reinv,  the  inverse  of  the  Reynolds  number 

par  (3)  alphar,  real  part  of  the  complex  wavenumber  in  direction 

pair  (4)  alphar,  imaginary  part  of  the  complex  wavenumber  in 
direction  a,- 

pair (5)  beta,  real  wavenumber  in  direction 

pair  (6)  omegam,  real  part  of  complex  temporal  frequency,  ujr 
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paa:(7)  omegai,  imaginary  part  of  complex  temporal  frequency, 

par  (8)  UsqdRT,  UURT^  =  7 

paLr(9)  coordinate,  i.e.,  azimuthal  amgle  4> 

6.1.9  Parameters 

Parameters  and  taisks  of  Linear  stability  analysis  to  be  performed  in  running  LISA 
are  assigned  by  the  user  in  the  input  file  Parameters.  The  first  example  in  Figure  6.3 
is  for  calculating  local  eigenvalues  and  eigenvectors  at  a  specified  location. 


1 

3  4 

2338.535 
4.2761798E-04 
2.1000000000E-01 
-3.656700000E-05 
O.OOOOOOOOOOE+01 
1.9340000000E-01 
O.OOOOOOOOOOE+01 
8.9600000000E+01 
O.OOOOOOOOOOE+01 
3 
1 

6  20  .004 

-999/END 


/task  : 

/eigenvalue  (set  to  alpha) 

/parameter  1  xi''{l)  (=re,  usually) 

/parameter  2  1/re 

/parameter  3  alphar 

/parameter  4  alphai 

/parameter  5  beta 

/parameter  6  omegar 

/parameter  7  omegai 

/parameter  8  U**2/RT 

/parameter  9  xi '“  { 3 } 

/task 

/kdim 

/parameter,  steps,  step  size 


Figure  6.3:  Example  Input  File  Parameters  for  Finding  Spatial  Eigenvalues  and 
Generating  a  Table  for  Varying  Wr  at  Station  s  =  s*/rl  =  2338.535  for  2D  Second- 
Mode  Disturbances. 

Here  itask=l  is  specified  as  a  first  task.  Parameters  3  and  4,  i.e.,  real  and 
imaginary  of  wavenumber  a  are  defined  as  unknown  eigenvalues.  The  streamwise 
location  is  specified  in  par(l)  in  dimensionless  units  (in  terms  of  boundary-layer 
length  scale  rl  =  \Jvs*  jUoo  at  the  reference  station  s*  JRf^  =  175).  The  inverse  of 
the  Reynolds  number  at  the  reference  station  is  given  in  par  (2).  The  eigenvalues 
peu:(3)  and  par  (4)  are  given  initial  guesses.  The  wavenumber  in  zizimuthal  direc¬ 
tion  pax (5)  is  set  to  zero  to  study  2D  disturbzinces.  The  dimensionless  frequency 
Ur  =  /*  rl/(27r  Uoo)  is  specified  as  par (6).  The  imaginary  part  w,  is  set  to  zero  for 
calculating  spatial  growth.  The  value  of  par  (8)  =  is  89.6  with  Moo  =  8-  The 
basic  state  is  axisymmetric  and  hence  par(9)=(^  is  set  to  zero.  After  this  task  is 
finished,  the  next  teisk  itask=3  generates  a  table  of  solutions  for  different  Ur-  The 
step  size  for  table  may  be  assigned  a  negative  value. 

The  second  example  (Figure  6.4)  first  performs  (itask=2)  at  the  initial  station  s  = 
s*/rl  =  1400.  Hence  eigenvalues  and  the  corresponding  eigenfunctions  are  calculated 
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2 

3  4 

1400.0 

4.276179o£:-04 

2.1000000000E-01 

-3.656700000E-05 

O.OOOOOOOOOOE+01 

1.9340000000E-01 

O.OOOOOOOOOOE+01 

8.960000000E+01 

O.OOOOOOOOOE+01 

3 

1 

1  78  30 

-999/END 


/task  : 

/eigenvalue  (set  to  alpha) 
/parameter  1  xi^{l}  (=re, usually) 
/parameter  2  1/re 
/parameter  3  alphar 
/parameter  4  alpha! 

/parameter  5  loeta 
/parameter  6  omegar 
/pareuneter  7  omegai 
/parameter  8  U**2/RT 
/parameter  9  xi { 3 ) 

/task 

/kdim 

/parameter,  steps,  step  size 


Figure  6.4:  Example  Input  File  Parameters  for  Finding  Spatial  Eigenvalues  and 
Generating  a  Table  for  Varying  =  s  at  the  Fixed  Frequency  of  Ur  =  0.1934  for  2D 
Second-Mode  Disturbances. 

at  this  station.  Eigenfunctions  axe  written  out  to  the  output  file  0. function.  Then 
a  series  of  solutions  axe  generated  (itask=3)  at  the  fixed  values  of  other  specified 
parameters  for  several  stations  by  varying  s.  The  file  0  .function  is  overwritten  each 
station  so  that  only  the  eigenfunctions  at  the  last  station  will  be  saved. 

6.1.10  Control 

An  example  input  file  Control  for  a  nonlinear  PSE  run  for  two-dimensional  distur¬ 
bances  is  given  in  Figure  6.5. 

Comments  provide  brief  explanations  of  entries  in  the  file.  There  eire  five  distinct 
modes  to  be  included  in  this  calculation.  The  initial  mode  is  the  fundamental  mode 
(1,0)  with  dimensionless  frequency  tpar.  Mode  control  mgr  for  this  mode  is  2,  the 
default  value.  This  mode  is  calculated  by  the  local  lineeir  stability  analysis  with 
itask=2  and  read  in  by  PSH.  Its  initial  amplitude  for  mjiximum  disturbance  (i.e., 
normv  is  set  to  3)  is  specified  as  2  x  10“^.  Three  other  modes  are  harmonics  of  this 
initial  mode,  to  be  created  by  the  PSH  through  subroutine  newmodes  if  the  nonlinear 
forcing  terms  are  greater  than  the  specified  thron  of  5  x  10“®.  The  phase- locking 
condition  will  be  used  for  these  harmonics,  as  mgr’s  axe  set  to  1.  The  fifth  (0,0) 
mode  is  meanflow  distortion,  mgr  for  this  mode  is  also  set  to  1. 
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Transition  in 

Hypersonic  Flow  over 

a  Blunt  Cone  at  Mach  8 

1  26 

1  20 

/Initial  and  final 

stations;  station  change;  step  size 

0 

250 

.1934 

0 

/inlt,  nterms  treated  implicitly  for  1,  explicitly  for  0 
/Iterations  (linear  <=  1) ,  max  step  divisions 
/tpar  (frequency),  zpar  (beta) 

1  0 
5e-9 

5e-10 

/itr,  izr  (reference  mode) . 

/Thresholds  on(rhs) /off  (eui^) 

14 

0  0 

2 

/Nximber  of  mode  control  lines 

/T- index,  Z-index,  include  this  mode. 

0  1 

-1 

/T- index,  Z-index, 

disregard  this  mode. 

0  2 

-1 

/T-index,  Z-index, 

disregard  this  mode. 

2  0 

1 

/T-index,  Z-index, 

include  with  phase  loclc. 

1  1 

-1 

/T-index,  Z-index, 

disregard  this  mode. 

1  2 

-1 

/T-index,  Z-index, 

disregard  this  mode. 

2  1 

-1 

/T-index,  Z-index, 

disregard  this  mode. 

2  2 

-1 

/T-index,  Z-index, 

disregard  this  mode. 

3  0 

1 

/T-index,  Z-index, 

include  with  phase  loclc. 

3  1 

-1 

/T-index,  Z-index, 

disregard  this  mode. 

3  2 

-1 

/T-index,  Z-index, 

disregard  this  mode. 

4  0 

1 

/T-index,  Z-index, 

include  with  phase  lock. 

4  1 

-1 

/T-index,  Z-index, 

disregard  this  mode. 

4  2 

-1 

/T-index,  Z-index, 

disregard  this  mode. 

1 

10  2 

2.e-4 

/Nvunber  of  initial  modes 

0./ T-index,  Z-index,  mgr,  amp,  phase 

Figure  6.5:  Example  Input  File  Control  for  the  Nonlinear  PSE  Run  for  Two- 
Dimensional  Second-Mode  Disturbance. 
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6.2  Hypersonic  Flow  over  a  Sharp  Cone  at  Angle 
of  Attack 

6.2.1  metric.ins 

The  disturbance  coordinate  system  is  nonorthogonal  and  almost  identical  to  the  al¬ 
gebraic  one  that  the  AFWAL  PNS  code  uses.  It  is  given  by: 

=  Rsin(f-ir)  =  -i?sin(f3) 
x^  =  Rcos^~n)  =  — i?cos({^)  ^  ' 

R  =  r,{i\e)+e 

where  is  the  distance  measured  along  the  body  axis,  is  the  circumferential  angle 
measured  in  radians  from  the  windward  meridian  in  the  direction,  is  the  local 
body  radius,  and  is  the  radial  distance  from  the  body  surface.  See  Figure  6.6.  The 
similarity  between  the  AFWAL  PNS  algebraic  grid  and  the  coordinate  system  used 
for  the  stability  analysis  is  intentional  and  enables  the  basic  state  to  be  interpolated 
with  greater  ease  and  accuracy. 

The  user  can  either  use  this  coordinate  system  or  a  local  Cartesian  coordinate 
system  oriented  tangent  to  the  body  surface  by  setting  the  10^,  or  lO’s  digit,  in 
blmodel  to  1  or  0,  respectively,  in  the  Definitions  file.  (This  digit  is  named  icoord 
in  the  metric .  ins  file.)  In  the  latter  case,  any  curvature  of  the  body  surface  will  be 
neglected. 


6.2.2  gridbs.ins 

Subroutine  gridbs  assigns  the  basic  state  grid  simply  by  reading  it  in  along  with 
the  basic  state  solution.  Auxiliary  routines  that  gridbs  calls  to  accomplish  this  are 
described  more  fully  in  Appendix  C.  Of  these,  three  axe  specific  to  the  basic  state 
flow  solver  used.  They  axe:  bsiois,  bpstda,  and  sc2G.e  . 

Subroutines  bsiois  and  bpstda  read  in  the  station  data  and  flow  profiles,  respec¬ 
tively,  from  FORTRAN  unformatted  PLOT3D  files.  These  files  contain  consecutive 
axial  planes  of  data  created  by  the  AFWAL  PNS  code. 

The  AFWAL  Parabolized  Navier-Stokes  (PNS)  code  has  been  developed  by  sev¬ 
eral  groups  since  1979  and  has  been  documented  by  Kaul  and  Chaussee  (1983);  Stal- 
naker,  Nicholson,  Hanline,  and  McGraw  (1986);  and  Rajendran  (1989).  It  uses  a 
space-marching  technique  and  body-fitted  coordinates  to  integrate  the  PNS  equa¬ 
tions  downstream  of  data  specified  in  an  initial  axial  plane  x^  =  xj. 

The  permissible  AFWAL  PNS  body-fitted  coordinates  axe  of  the  form: 

x^  =  x*(/) 
x^  =  x\J,K) 
x^  =  x%J,K) 


51 


(6.2) 

(6.3) 

(6.4) 


where  is  the  body  axis  (see  Figure  6.6).  This  means  that  the  solution  is  marched 


Figure  6.6:  Coordinate  System  Used  in  the  Stability  Analyses. 

in  planes  perpendicular  to  the  body  axis.  Hence,  the  body-fitted  coordinates  are 
nonorthogonal  whenever  the  slope  of  the  body  in  the  axial  direction  is  nonzero.  This 
is  certainly  the  case  for  any  cone  of  nonzero  half-angle. 

Virtually  the  same  coordinate  system  (Section  6.2.1)  is  used  for  the  stability  aneil- 
yses  to  simplify  the  problem  of  interpolating  the  beisic  state  solution  onto  the  dis¬ 
turbance  grid.  The  interpolation  scheme  does,  however,  place  restrictions  on  the 
transformation  between  the  basic  state  coordinates  {K,  I,  J)  and  the  disturbance  co¬ 
ordinates  These  restrictions  axe  discussed  in  Appendix  C  . 

The  PLOT3D  files  must  be  created  according  to  the  formats  used  by  the  following 
routines  to  read  them. 

bsiois 

Subroutine  bsiois  transforms  the  Cartesian  grid  points  in  the  body  surface  to  cylin¬ 
drical  coordinates  and  stores  the  required  data  in  the  matrix  stndat  . 

The  relevant  open,  read  and  station  data  assignment  statements  are: 

one  =  1. 

pi  =  acosC-one) 

C  Open  the  PL0T3D  grid  file  — 
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open(uiiit«ub8f (1) ,fil«>b8fn(l),acce8s«’sequential’ ,8tatus»’old’ , 

I  form* ’ unf ormatt«d ’ ) 

r«ad(ub8f (1))  idim, jdim.kdim 

C  Open  the  PL0T3D  solution  file  . . . 

open(unit*ub8f (2) ,file*b8fn(2) .access* ’sequential’ ,status»’old’ , 

I  form* ’unformatted’) 

read(ubsf (2))  idl.id2.id3 

if  (idl.ne.idiffl  .or.  id2.ne.jdiffi  .or.  id3.ne.kdiffl)  then 
write(6.98)  idim.jdim.kdim.idl.id2.id3 
98  format(//,’  ’ ,72(’*’) .//.3X. ’The  dimensions  of  the  basic  ’. 

I  ’state  grid  and  solution  files  eo-e  not  the  same.’.//. 

I  3X.’GRID:’.5X.I3.’  X  ’.13.’  x  ’ .I3./.3X. ’SOLUTIOH: ’ . 

I  IX. 13.’  X  ’.13.’  X  ’.13.//.’  ’.72(’*’)) 

stop 
end  if 
C 

C  Free-stream  Re  assumed  to  be  based  on  free-stream  speed  of  sound  and 
C  NOT  the  free-stream  flow  speed, 
c 

read(ubsf (2))  fsmach. attack. fsre.fstemp 


istn*0 

nk*0 

do  300  k*l.kdim 

readCubsf (1) ,iostat*ie(l))  ((xb  .i*l,idim) . j*l, jdim)  , 

I  ((yb(i.j) .i*l,idim) .j*l,jdim) , 

I  ((zbCi.j) .i*l,idim) .j*l,jdim) 

C 

C  Don’t  repeatedly  store  data  for  the  same  station  ... 

C 

if  (k.gt.l.and.xb.eq.xbo)  go  to  400 

if  (ieCl) .eq.0.and.ie(2) .eq.O)  then 

xbo  *  xb 
nk  *  nk+1 

if  (k.eq.l)  then 
C 

C  Teirget  value  for  XI3BSC  computed  from  YBCl.l)  and  ZB (1.1).  The 
C  branch  used  for  the  arctangent  will  be  chosen  so  that  XI3(K*1.J*1) 
C  will  be  closest  to  this  tairget  value. 


53 


c 


903 


write(6,903)  xb,yb(l,l) .zbCl.l) 

formate/,’  ’Converting  PL0T3D  grid  file  from  Cartesian’ 
,’  coordinates  to  cylindrical  coordinates.’,/,’  ’, 
’Hust  choose  branch  of  arctangent.  Enter  the  ’, 
’circumferential  angle  for  point  NEAR’,/,’  ’ , ’ (not  ’, 
’  necessarily  coincident  with)  the  point:’,/,’  ’, 

’XB  «  ’,E12.6,’  YB(1,1)  =  ’,E12.6,’  ZB(1,1)  =  ’, 
E12.6,’  =>  ’,$) 
read(5,*)  xi31st(0) 
xi31st(0)  *  pi 
w  s  one 
else 


C 

C  For  k>l,  choose  the  branch  of  the  arctangent  so  that  the  computed 
C  value  of  xi3(k,l)  is  near  xi3(k-l,l). 

C 

xi31st(0)  =  xi31st(l) 

V  =  0.5 


do  100  j=l,jdim 
istn  =  istn+1 
xilbsc  =  xb 

xi3bsc  ®  atan2(-yb(l,j) ,-2b(l, j)) 

C 

C  Choose  the  branch  of  the  arctangent  so  that  xi3bsc  is  near  the 
C  average  value  of  xi3(k-l,j)  (stored  in  xi31st(j))  and  xi3(k,j-l) 

C  (stored  in  xi31st(j-l)) . 

C 

branch  =  nint((w*(xi31st(j)+xi31st(j-l))-xi3bsc) 

I  /(2.*pi)) 

xi3bsc  =  xi3bsc  +  2*branch*pi 

if  (k.eq.l)  write(6,*)  ’xi3bsc=’ ,xi3bsc, ’  branch=’ , branch 


stndat(  0,istn) 
stndat(  l,istn) 
stndat(  2, istn) 
stndat(  3, istn) 
stndat(  4, istn) 


1+idim* ( (k- 1 ) ♦ j dim+ ( j - 1 ) ) 

xilbsc 

idim 

xi3bsc 

k 


C  Local  body  radius  . . . 


stndat(  5, istn) 


sqrt (yb ( 1 , j ) ♦♦2+zb ( 1 , j ) **2) 


C  Local  shock  stand-off  distance  measured  from  body  axis  . . . 
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stndatC  6,istn)  *  8qrt(yb(idiiB, j)**2+zb(idim, j)**2) 

xiSlstCj)  s  ziSbsc 
100  continue 

else 

C  Write  diagnostics  for  read  errors  . . 

end  if 

300  continue 

400  close(ubsf (1)) 
closeCubsf (2)) 

where 

xilbsc  = 
xiSbsc  = 
yb  = 
zb  = 

idim  =  number  of  basic  state  grid  points  in  ^ 

j  dim  =  number  of  basic  state  grid  points  in  ^ 

kdim  =  number  of  basic  state  grid  points  in 

Typically,  the  AFWAL  PNS  code  will  store  the  flow  field  in  the  range  tt  <  < 

27r  when  the  flow  and  body  axe  symmetric  about  the  =  z  axis.^  For  this  half¬ 
interval,  the  user  should  enter  a  value  of  neeir  3  w  tt  to  select  the  right  branch  of 
the  arctangent  needed  to  convert  the  PLOT3D  Cartesian  coordinates  to  cylindrical 
coordinates. 

Subroutine  bsiois  also  sets  the  sc2des  for  the  basic  state  solution.  The  PLOT3D 
bsiois  sets  the  reference  length  to  1”  for  and  It  then  converts  this  reference 

length  to  units  of  meters  before  storing  it.  The  scale  for  the  azimuthal  angle 

measured  aroimd  the  body  axis,  is  simply  set  to  1.  Also,  bsiois  assumes  that 
the  reference  temperature  for  the  basic  state  cedculation  is  assigned  to  the  variable 
normally  reserved  for  time  on  the  second  line  in  the  PLOT3D  file: 

^We  did  not  realize  this  at  the  beginning  of  our  anedyses  with  the  PNS  code.  In  retrospect, 
it  would  have  been  better  to  define  =  0  along  the  windward  meridizui  so  that  the  computed 
half-interval  would  be  0  <  <  z-.  This  could  be  dianged  by  modifying  the  file  metric .  ins  . 
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C  Open  the  PL0T3D  solution  file  — 

open(unit*ubsf (2) ,f ile«b8fn(2) .access^’ sequential’ .statuss’old’ , 

I  f oxm« ’ unf ornatted ’ ) 

readCubsf (2))  idim. jdis.kdim 
C 

C  Free-stream  Re  assumed  to  be  based  on  free-stream  speed  of  sound  and 
C  NOT  the  free-stream  flos  speed, 
c 

r ead (ubsf ( 2 ) )  f smach , att  ack , f sre . f st  emp 

Although  this  temperature  is  assumed  to  be  measured  in  degrees  Rankine  (re:  AFWAL 
PNS),  subroutine  bsiois  stores  it  in  Kelvin.  Finally,  the  Reynolds  number  which  is 
read  in  is  based  on  the  free-stream  speed  of  sound,  but  subroutine  bsiois  stores  the 
Reynolds  number  based  on  the  free-stream  velocity.  All  of  these  conventions  can  eas¬ 
ily  be  changed  by  modifying  the  corresponding  assignment  statements  in  subroutine 
bsiois  . 


bpstda 

Subroutine  bpstda  transforms  the  Cartesian  grid  to  cylindrical  coordinates  and  the 
Cartesian  velocity  components  to  the  contravariant  components  in  the  stability  co¬ 
ordinates.  It  also  computes  the  temperature  from  the  total  energy  per  unit  volume 
that  is  read  in  and,  finally,  writes  out  all  the  transformed  data  to  the  direct  access 
file  station  by  station. 

The  relevant  open  and  read/write  statements  for  the  PLOT3D  files  are:^ 
parameter  (ibsl=501 ,ibs2-19,nbsvar-5) 


gamma  -  realvClS) 

gnwfil  s  gSJnOBlSl^OZIO 
gr^rnml  s  g3]0lU3l^^IDnLl 


open(unit=ubsf (1) ,file=bsfn(l) ,access=*sequential’ ,status=’old’ , 
I  f orm= ’ unformatted ’ ) 

open(unit=ubsf (2) ,file-bsfn(2) .access* ’sequential’ ,status=’old’ , 
I  form®’unformatted’ ) 


readCubsf (1) ,iostat*ie(l))  idim, jdim.kdim 
readCubsf (2) ,iostat=ieC2))  idl,id2,id3 

if  Cidl.gt.ibsl  .or.  id2.gt.ibs2)  then 
writeC6,99)  idl,id2,ibsl,ibs2 

99  formate//,’  ’ ,72C’*’) ,//,3X, ’The  dimensions  of  the  basic  ’, 
^The  reed  variables  in  the  PLOT3D  files  must  be  single  precision. 
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o  o 


I 

I 

I 


’state  grid  and  solution  files  are  too  big.’,//, 

3X. 'FILES:  ’,5X.I3,’  x  ’ .13,/, 3X. 'MEMORY: ’ , 

IX. 13.’  X  ’.13.//,’  ’,72(’*’)) 

stop 
end  if 

readCubsf  (2)  ,iostat>ie(2))  fsmacb., attack, fare, fsta&p 
istn»0 

do  300  k=l,kdim 

readCubsf (1) ,iostat»ie(l))  ((xb  ,i*l,idim) , j»l, jdim) , 

I  ((yb(i,j),i=l,idia),j*l,jdiin), 

I  ((zb(i,j),i=l,idim) ,j»l,jdim) 

readCubsf  (2)  ,iostat»ieC2))  CCCqbCi.j  ,nv)  ,i=l,idiiii) ,  j=l,  jdim)  , 
I  nv®l,nbsvar) 


if  CieCl) .eq.0.and.ieC2) .eq.O)  then 

do  120  j®l,jdim 
istn  =  istn+1 

key  3  nintCstndatCO.istn)) 

Rotate  into  cylindrical  coordinates  Ccontravariant  components) . 

C  Assume  that  the  radial  lines  are  STRAIGHT! 

C  This  step  is  essentially  the  stability  anailysis  grid  generation. 

C  Also,  scale  the  velocity  components  with  the  free-stream  flou  speed 

C  instead  of  the  free-stream  speed  of  sound. 

C 

kbamiaxCminCk.kdim)  ,2) 

dlrdl®  CstndatCS, Ckb-l)*jdim+j)-stndatC5, Ckb-2)*jdim+j)) 

I  / Cstndat  C 1 , Ckb- 1 ) * j  dim+j ) -stndat  C 1 , Ckb-2) * j  dim+ j ) ) 

dlrdS  =  0 

do  110  i®l,idiffl 

C  eta  *  atan2C-ybCi,j),-zbCi,j)) 

rb  =  sqrtCybCi,j)*ybCi,j)+zbCi,j)*zbCi,j)) 
bvelCl)  *  qbCi.j ,2)/qbCi.j,l)/fsmach 
bvelC2)  =  qbCi,j,3)/qbCi,j,l)/fsmach 
bvelC3)  =  qbCi.j ,4)/qbCi,j,l)/fsmach 
qbCi,j,2)  ®  Cggaiinl/qbCi,j,l))*CqbCi,j,5)-0.5d0/qbCi,j  ,1) 
I  *CqbCi,j,2)**2+qbCi,j,3)**2+qbCi,j ,4)**2)) 

C  Contravariant  velocity  components  in  cylindriced  coordinates  . . . 
qbCi.j ,3)  *  bvelCl) 
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qb(i,j,4)  ■  (yb(i,j)*bvel(2)+zb(i,j)*bvel(3))/rb 
qb(i,j,5)  ■  (zb(i,j)*bvel(2)-yb(i, j)*bvel(3))/rb/rb 

C  Contravariant  velocity  components  in  cylindrical  coordinate  system 
C  vith  zi"2  »  0  at  surface  . . . 

qb(i,j,4)=  -dlrdl*qb(i,j,3)  +  qb(i,j,4)  -  dlrd3*qb(i, j ,5) 

write(uda(l) ,rec*key+i-l)  rb, (qb(i, j ,nv) ,nv»l,5) 

110  continue 

120  continue 
end  if 

300  continue 

6.2.3  bstate.ins 

Subroutine  bstate  reads  in  and  interpolates  the  basic  state  profiles  from  the  direct 
access  file  created  by  subroutine  gridbs.  The  method  by  which  it  does  this  is  inde¬ 
pendent  of  the  input  and  output  routine  required  to  read  in  the  baisic  state,  and  so 
need  not  be  changed  for  different  basic  state  flow  solvers.  It  is  described  more  fully 
in  Appendix  C. 

6.2.4  btrans.ins 

This  file  specifies  the  relationship  between  the  basic  state  variables  that  axe  used 
in  the  stability  equations  and  the  basic  state  variables  that  axe  used  in  subroutine 
bstate  .  As  described  in  Section  2.2.4,  this  relationship  may  be  selected  by  specifying 
density,  temperature,  and  one  of  three  sets  of  velocity  components.  These  three  sets 
of  velocity  components  are; 

1.  Cartesian  components  in  an  inerti2d  system. 

2.  Physical  contravariant  components  in  the  disturbance  coordinates. 

3.  Contravariant  components  in  the  disturbance  coordinates. 

The  user  specifies  the  set  of  velocity  components  for  btrans  by  setting  the  integer 
variable  ibvel  =  0,1,  or  2,  respectively. 

Although  the  AFWAL  PNS  solution  and  the  PLOT3D  files  that  it  creates-  use 
Cartesian  velocity  components,  the  velocity  vector  in  btrans  is  specified  in  te;ms 
of  the  contravariant  components  (ibvel  =  2),  because  the  AFWAL  PNS  code  runs 
axisymmetric  calculations  fully  three-dimensional  with  a  coarse^  point  distribution 

^Assuming  that  the  reference  coordinates  used  for  the  PNS  solution  are  cylindrical  -  an  option 
that  the  AFWAL  PNS  user  can  select. 
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in  the  circumferential  direction.  The  Cartesiein  velocity  components,  since  they  vaxy 
sinusoidally  around  the  circumference  of  an  axisymmetric  body  at  zero  tingle  of  at¬ 
tack,  would  not  be  interpolated  very  accurately  on  a  coarse  grid.  The  contravaxiant 
components,  which  are  invariant  in  the  circumferential  direction,  can  be  interpolated 
accurately. 

Two  different  options  are  also  implemented  to  control  some  of  the  terms  in  the 
stability  equations.  The  contravaxiant  velocity  component  and  all  and  deriva¬ 
tives  of  the  basic  state  can  be  set  identically  to  zero  by  setting  the  10°,  or  I’s  place, 
in  blmodel  to  zero  in  the  Definitions  file,  (blmodel  =  0 . )  Alternatively,  all  of  these 
terms  are  retained  in  the  analysis  for  a  nonzero  value  of  this  same  digit  in  blmodel . 

Subroutine  btrans  also  sets  the  the  thermophysical  properties  for  the  disturbance 
analyses.  This  is  done  through  calls  to  subroutine  therm,  described  in  Appendix  A  . 
The  user  can  select  which  constitutive  equations  therm  must  use  for  the  properties  by 
setting  the  second  constant,  Air-Prop,  in  the  Definitions  file  to  a  value  as  described 
in  Appendix  A  . 

6.2.5  boundary.ins 

The  boundary  conditions  which  axe  implemented  in  boundaury .  ins  axe  homogeneous. 
No  options  enable  the  use  of  asymptotic  boundary  conditions  in  the  fax-field,  but 
this  could  be  changed.  Also,  the  disturbance  is  also  assumed  to  satisfy  the  no-slip 
and  kinematic  boundary  conditions  at  the  body  surface  and  the  disturbance  in  the 
derivative  of  the  temperature  with  respect  to  is  set  identically  to  zero. 

6.2.6  pointd.ins 

The  radial  point  distribution  can  be  made  to  vaxy  with  by  modifying  the  grid 
paxameters  in  pointd.ins  .  The  current  setting  makes  chy(jx),  the  point  furthest 
from  the  surface,  vary  according  to: 

btrf  =  strfp(2)  ♦  sqrt(pau:(l)  *  pEu:(2))  (6.5) 

i.e.,  like  the  square  root  of  the  ratio  of  where  is  the  value  of  at  the  station 
where  the  disturbance  reference  length  =  ^o*/ 

6.2.7  post  pr. ins 

No  postpr.ins  file  was  written  for  this  application.  However,  the  corresponding  file 
for  the  blunt  cone  application  (Section  6.1.7)  could  be  modified  and  used  here  as  well. 
Data  for  the  figures  shown  in  the  final  report  were  taken  directly  from  the  0. table 
file  and  0 .  summary  file  created  by  the  local  stability  code  and  PSE  code,  respectively. 
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6.2.8  Definitions 


The  Definitions  file  is  virtually  identical  to  the  one  described  for  the  blunt  cone 
application  (Section  6.1.8).  The  differences  are  in  the  numerical  values  of  the  refer¬ 
ence  data  that  axe  input  and  the  controls  for  the  boundary  conditions.  There  are  no 
controls  for  the  sharp  cone  boundary  conditions  (see  Section  6.2.5).  The  print  control 
constant  npprt  is  also  not  present  because  a  postpr .  ins  file  was  not  written  for  the 
sharp  cone  application.  Finally,  the  meanings  of  the  parameters  are  also  different 
because  the  coordinate  systems  are  different.  (See  Sections  6.1.1  and  6.2.1,  respec¬ 
tively.)  While  the  specific  meaning  of  the  parameters  emd  some  constants  is  different 
for  different  coordinate  systems,  their  general  definitions  are  the  same  if  interpreted 
in  terms  of  aribtrary  (^^,^^,^^). 

6.2.9  Parameters  and  Control 

The  Parameters  and  Control  files  can  be  specified  for  the  sharp  cone  application 
just  as  they  were  for  the  blunt  cone  application.  See  Sections  6.1.9  and  6.1.10,  re¬ 
spectively. 
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Appendix  A 

Constitutive  Equations 


The  constitutive  relations  close  the  system  of  flow  equations  by  specifying  how  the 
properties  of  the  gas  depend  on  its  state.  A  variety  of  relations  for  each  required 
property  of  air  are  implemented  in  subroutine  therm,  therm  computes  the  dimen¬ 
sional  properties,  scales  them,  amd  then  returns  dimensionless  data.  There  is  no  need 
to  compute  and  specify  a  long  list  of  dimensionless  constants.  Instead,  the  input  of  a 
short  list  of  information  is  sufficient  to  determine  all  the  scales. 

The  user  selects  the  relations  and  nondimensionalization  by  assigning  a  value  to 
a  seven  digit  integer  model  number  named  mdlair.  mdlair  is  passed  as  the  first 
argument  to  the  subroutine  therm.  Each  digit  of  mdlair’s  base  10  representation  is 
itself  a  model  number  for  a  particular  constitutive  relation.  Hence,  there  may  be  up 
to  10  different  models  for  each  property. 

The  lowest  order  digit,  10°,  or  I’s  place,  in  mdlair  specifies  how  to  compute  the  scales 
used  to  nondimensionalize  the  properties.  Currently,  any  of  four  different  choices  may 
be  selected.  They  are  all  biased  towards  compressible  flow  applications  because  they 
all  require  that  the  reference  temperature  and  either  the  Mach  number  or  also  be 
input.  With  this  additional  data,  and  assuming  that  in  its  reference  state  air  behaves 
as  an  ideal  gas,  therm  can  compute  the  reference  speed  of  sound,  reference  velocity 
and  reference  viscosity  /x^.  How  therm  proceeds  to  compute  the  remaining  scales 
depends  on  which  of  the  four  options  the  user  selects  from  Table  A.l. 

The  first  two  options  in  Table  A.l  are  available  for  those  situations  when  most,  but 
not  all,  of  the  dimensionless  parameters  are  known.  These  are  important  caises  be¬ 
cause  all  the  properties  that  therm  computes,  except  for  the  density,  vary  only  with 
temperature.  Subroutine  therm  can  thus  compute  the  (dimensional)  reference  values 
of  all  the  properties  except  for  the  density  knowing  only  the  reference  temperature. 
Using  these  reference  values  as  scales,  therm  cjin  subsequently  nondimensionalize  all 
the  constants  appearing  in  all  the  constitutive  equations  except  for  the  equation  of 
state.  The  equation  of  state  for  an  ideal  gais,  on  the  other  hand,  can  easily  be  coded 
in  dimensionless  form  and  depends  on  only  one  dimensionless  constant:  U^I{R*T^). 
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10°  digit  Selection 


0 

Input  the  reference  pressure  (atm)  and  the  inverse  of  the 
Reynolds  number.  therm  computes  the  reference  density 
(kglm^)  from  the  ideal  gas  law  for  air.  It  then  determines  the 
reference  length  (m)  from  the  inverse  of  the  Reynolds  number. 

1 

Input  the  reference  length  (m)  «ind  the  inverse  of  the  Reynolds 
number,  therm  computes  the  reference  density  (kg/m^)  from 
the  inverse  of  the  Reynolds  number.  It  then  determines  the 
pressure  (atm)  from  the  ideal  gas  law  for  air. 

2 

Input  the  reference  pressure  {atm)  and  the  reference  length  (m). 
therm  computes  the  reference  density  {kg/m^)  from  the  ideal 
gas  law  for  air.  It  then  computes  the  inverse  of  the  Reynolds 
number. 

3 

Input  the  reference  pressure  (atm)  and  a  length  c  (m).  therm 
computes  the  reference  density  {kgfm^)  from  the  ideal  gas  law 
for  air.  It  then  computes  the  inverse  of  the  Reynolds  number 
Rcc  based  on  c.  Finally,  the  reference  length  is  set  to  cfy/Rcc 
and  a  new  Reynolds  number  =  y/Rec  based  on  this  reference 
length  is  computed. 

Table  A.l:  Scales  selection 
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10*  digit 

Selection 

0 

^  13  input 

1 

Mach  number  is  input. 

Table  A. 2:  Parameter  selection 

Hence,  for  fixed  values  of  this  constant^  and  for  fixed  reference  temperature,  «dl  the 
dimensionless  constitutive  equations  are  independent  of  the  reference  pressure  or  ref¬ 
erence  length  that  the  user  inputs  and  consequently,  so  is  the  dimensionless  solution 
for  the  flow. 

The  options  that  set  the  10°  digit  =  0  or  1  are  typically  used  when  studying  the 
stability  of  the  flat  plate  compressible  boundary  layer  with  similar  velocity  profiles. 
The  options  that  set  the  10°  digit  =  2  or  3  are  used  when  the  user  does  not  know 
the  Reynolds  number  and  wants  therm  to  compute  the  value  by  using  the  chosen 
constitutive  equations. 

Continuing  with  the  other  digits  in  mdlair,  the  10^  digit,  or  lO’s  place,  indicates 

whether  or  the  Mach  number  is  input.  See  Table  A.2. 

^  -*00 

The  10^  digit,  or  lOO’s  place,  indicates  the  equation  of  state.  The  equation  of  state 
specifies  the  density  as  a  function  of  pressure  and  temperature.  The  density  is  nondi- 
mensionalized  by  the  reference  density,  the  temperature  is  nondimensionalized  by  the 
reference  temperature,  and  the  pressure  is  nondimensionalized  by  the  reference  den¬ 
sity  multiplied  by  the  square  of  the  reference  velocity.  Table  A. 3  gives  the  different 
models  that  the  user  can  select  from. 

The  10^  digit,  or  lOOO’s  place,  indicates  the  constitutive  equation  for  the  specific  heat 
at  constant  pressure.  The  specific  heat  is  nondimensionalized  by  the  gas  constant. 
The  three  models  which  have  been  implemented  axe  given  in  Table  A.4. 

The  10^  digit,  or  10,  OOO’s  place,  indicates  the  model  for  the  first  coefficient  of  viscosity. 
Table  A. 6  identifies  the  four  different  relations  which  have  been  implemented.  The 
dynamic  viscosity  is  nondimensionalized  with  its  value  at  the  reference  temperature. 

The  10®  digit,  or  100,000’s  place,  indicates  the  model  for  the  second  coefficient  of 
viscosity.  Currently,  only  two  models  have  been  implemented.  See  Table  A. 7.  It 
is  nondimensionalized  by  the  value  of  the  first  viscosity  coefficient  evaluated  at  the 
reference  temperature. 

^  which  is  either  input  or  directly  calculable  knowing  the  input  Mach  number  and  reference 
temperature 
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10^  digit 

Selection 

0 

p*  =  incompressible  flow. 

1 

p*  =  ideal  gas.  The  gas  constant  i?*  used  for  air  is  287 

Jlikg.K). 

Table  A. 3:  Equation  of  state  selection 


10^  digit 

Selection 

0 

c*  =  where  7  =  1.4  is  a  constant.  (Calorically  perfect  gas.) 

1 

c!  =  as  +  air;,  +  air;2  +  a5r;3  + 

the  dimensional  constants  having  been  determined  by  Fabio 
Bertolotti  using  a  least  squares  fit  to  a  set  of  experimental  data 
in  the  temperature  range  lOOA’  <  T*  <  1600/^.  The  dimen¬ 
sional  constants  for  T*  in  K  and  Cp  in  J f  {kg  •  K)  axe  listed  in 
Table  A.5. 

Table  A.4:  Specific  heat  constitutive  relations 


ID 

b*  ip*,N  ■  s/m^) 

cUK^WIim-K)) 

ID 

ID 

! 

1 

1 

■iihiiikimriBHi 

ID 

-6.9301497  X  lO"” 

-6.8469791  x  lO"® 

ID 

Bli  dALW^dlliyiB 

3.3270833  x  10"^^ 

4 

1.9108582  X  10-^° 

j 

i 

I 

Table  A. 5:  Least  squares  polynomieil  coefficients  for  thermophysical  data 
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10^  digit 

Selection 

0 

/i*  =  constant  dynaimic  viscosity. 

1 

//*  =  linear  variation. 

2 

fi*  =  fi*  (r*/r*i^I/r»)‘  Sutherland’s  law,  where 

c;  =  110.4  K  and  fi;  =  1.719x10-®  •  s/m^  at  T;  =  273.1K. 

3 

^- =  6;  +  4;^  + 

the  dimensional  constants  having  been  determined  by  Fabio 
Bertolotti  using  a  least  squares  fit  to  a  set  of  experimental  data 
in  the  temperature  range  lOOA”  <  T*  <  IGOOA".  The  dimen¬ 
sional  constants  for  T*  in  K  and  fi  in  N  •  s/m^  are  listed  in 
Table  A.5. 

Table  A. 6:  Constitutive  relations  for  the  first  coefficient  of  viscosity 


Table  A. 7:  Constitutive  relations  for  the  second  coefficient  of  viscosity 
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10®  digit 

Selection 

0 

K*  =  =:  constant  thermal  conductivity. 

^oo^poo 

1 

K*  =  =:  constant  Premdtl  number. 

2 

K*  =  «r  Sutherland’s  law,  where 

=  194.4  a: 'and  <  =  \mxlQ-^WI{m  K)  bXT;  =  2n.\K. 

3 

K’  =  ci  +  ciTi  +  +  4T;i  + 

the  dimensional  constants  having  been  determined  by  Fabio 
Bertolotti  using  a  least  squares  fit  to  a  set  of  experimental  data 
in  the  temperature  range  lOOAT  <  T*  <  IGOOA.  The  dimen¬ 
sional  constants  for  T*  in  K  amd  k  in  Wj{m  •  K)  are  listed  in 
Table  A. 5. 

Table  A.8:  Constitutive  relations  for  the  thermal  conductivity 

Finally,  the  10®  digit,  or  1,000,000’s  place,  indicates  the  model  for  the  thermal  con¬ 
ductivity.  The  three  choices  are  given  in  Table  A.8.  The  thermal  conductivity  is 
nondimensionalized  by  the  product  of  the  gas  constant  and  the  value  of  the  first 
coefficient  of  viscosity  evaluated  at  the  reference  temperature. 
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Appendix  B 

Asymptotic  Boundary  Conditions 


The  use  of  asymptotic  boundary  conditions  in  solving  the  stability  equations  (PSH 
and  LSH)  allows  the  farfield  boundary  to  be  located  closer  to  the  wall.  Hence  fewer 
grid  points  are  required  for  adequate  resolution.  It  also  eliminates  oscillations  in  the 
phaise  of  the  shape  function  at  the  farheld  boundary  observed  in  solutions  obtained 
with  homogeneous  boundary  conditions  at  a  tnmcated  normal  distance.  The  bound¬ 
ary  conditions  are  based  on  the  local  decay  rate  in  normed  direction  of  the  inviscid 
solutions  (Mack  1984).  The  mean  flow  velocity  is  assumed  to  be  uniform  locally, 
which  is  a  reasonable  assumption  provided  the  boundary  is  located  away  from  the 
shock  and  about  1.5  time  the  boundary-layer  thickness  from  the  wall. 

In  body-intrinsic  coordinates 

=  s;  arc-length  along  the  body 
=  Tj]  normal  distance  from  the  body  surface 
=  <f>-,  azimuthal  angle 

Wave  numbers  are  defined  as 

ds 
dz 

d<f> 


m 

dt' 

The  normal-mode  solution  for  pressure  disturbance  in  inviscid  uniform  freestream  is 
given  in  Mack  (1984).  It  behaves  like 


a 

and  the  frequency  is  defined  as 
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where  Mg  is  the  relative  (complex)  Mach  number  in  the  freestream  given  by 

^  _  {cxU,  +  ^W,-u)M, 

Here  U  and  W  be  mean-flow  velocity  in  streamwise  x  and  z  directions,  respectively. 
The  subscript  e  denotes  freestream  values.  Other  disturbance  variables,  T,  u,  v  and 
w  should  decay  (in  y)  with  the  same  exponential  rate. 

The  above  equations  are  are  derived  by  Mack  in  Cartesian  coordinates.  A  Carte¬ 
sian  coordinate  may  be  constructed  locally  at  the  cone  surface  so  that  (r,  y,z  =  0) 
directions  align  with  (s,  t}^  <j))  of  the  axisymmetric  body  intrinsic  coordinates.  Here, 
6z  =  rS^  is  the  local  rectilinear  direction  and  r  is  the  r2wlius  from  the  cone  axis  to 
the  boundary  point.  The  wave  number  (integer  number  of  waves)  to  be  used  with 
<t>  can  be  related  to  the  rectilinear  wavenumber  ^  by  the  following  equation. 

dd  dz 

~  Tz’d^ 

=  {ri  +  r)co${$c)}.0 

If  =  d4>jdt,  then  the  rectilinear  velocity  =  {rj,  -|-  r]cosOc}w,i>  decays  exponen¬ 
tially  in  T)  with  the  same  rate  as  w.  Hence,  we  can  derive  the  following  differential 
equation  for  /  =  {T,  to  be  applied  as  a  boundau-y  condition. 

1^  +  0/  =  0 

OT] 

where  /  is  a  4  x  4  identity  matrix  and  is  1  for  f*  =  cind  is  zero  otherwise. 
The  equation  for  /  may  be  transformed  into  that  for  Q  other  desired  variables  in  the 
insert  file  boundauiy.  ins.  It  is  discretized  using  hermitian  finite-difference  formulae 
and  applied  at  the  bound^y  truncated  at  T^e  ~  1-5  —  where  is  the  boundary- 
layer  thickness.  For  the  hypersonic  boundary  layer  over  the  cone,  gradients  of  basic- 
state  variables  with  respect  to  rj  are  small  at  this  location.  The  error  in  the  boundary 
conditions  caused  by  the  nonuniformity  of  the  freestream  in  rj  should  be  small  as  it 
is  scaled  by  a  factor  of  0{ 
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Appendix  C 

Bstate.ins  for  Hypersonic  Cone 
Applications 


The  bstate  file  which  reads  the  FORTRAN  unformatted  PL0T3D  planes  files  or  the 
TLNS  solution  of  Esfahanian  (1991)  is  complicated  because  it  interpolates  a  basic 
state  solution  defined  numerically  on  a  grid.  It  currently  assumes  that  the  grid  is 
structured,  but  the  number  of  grid  points  along  the  coordinate  lines  passing  through 
the  body  surface  is  allowed  to  vary  from  station  to  station  (as  it  does  in  the  Kaups  & 
Cebeci  (1977)  boundary  layer  code).  Moreover,  the  bstate  file  currently  uses  tensor 
product  interpolation  (i.e.,  compositions  of  one-dimensional  interpolants)  and  works 
only  when  the  transformation  between  the  basic  state  coordinates  {K,  7,  J)  and  the 
disturbance  coordinates  (^^,^^,^^)  is  of  the  form: 


<■  =  cm  (c.i) 

e  =  (C.2) 

e  =  ew  (C.3) 

Some  minor  modifications  would  enable  the  treinsformation  to  be  of  the  form: 

e  =  e{K,j)  (C.4) 

e  =  e{K,i,j)  (C.5) 

e  =  e(K,j)  (C.6) 


These  modifications  have  not  yet  been  implemented.  More  extensive  modifications 
would  be  required  to  enable  the  transformation  to  be  arbitrary. 

The  restriction  that  and  be  independent  of  I  means  that  the  coordinate  lines 
must  have  the  same  shape  and  direction  as  the  I  coordinate  lines.  However,  each 
coordinate  line  is  labeled  with  the  values  of  (^^,^^)  which  are  constcmt  along  it 
instead  of  the  values  (i^,  J)  attached  to  the  corresponding  I  coordinate  line.  Hence,  if 
the  I  coordinate  lines  are  all  straight  and  perpendicular  to  the  body  surface,  so  must 
be  the  coordinate  lines.  If  the  I  coordinate  lines  are  all  straight  and  oblique  to  the 
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surface,  so  must  be  the  coordinate  lines.  Etc.  This  restriction  is  not  severe  because 
the  stability  codes  can  use  the  same  nonorthogonal  coordinates  that  modern  basic 
state  codes  use.  Also,  to  maie  the  boundary  layer  approximation  or  the  thin  layer 
approximation,  the  basic  state  viscous  flow  solvers  use  grids  with  the  I  coordinate 
lines  orthogonal  or  nearly  orthogonal  to  the  body  surface.  Hence,  no  adverse  effects 
on  the  accuracy  of  the  PSE  approximation  should  be  expected. 

The  interpolation  that  bstate  itself  is  responsible  for  is  primarily  along  the  coordi¬ 
nate  lines.  Another  subroutine,  mdintr,  determines  the  stations  and  interpolants  for 
interpolating  the  basic  state  in  and  (See  below.)  First,  bstate  calls  mdintr  . 
Next,  it  distributes  the  points  chy(j)  along  each  coordinate  line  specified  by 
mdintr  and  computes  the  corresponding  vatlues  of  the  basic  state  using  Lagrange 
polynomials  as  interpolants.  It  then  differentiates  the  Lagrange  interpolants  with 
respect  to  to  get  the  derivatives  of  the  basic  state  at  the  same  stations  and 
points.  Finally,  knowing  the  data  at  the  proper  points  but  at  surrounding  and 

stations,  bstate  interpolates  the  profiles  and  their  derivatives  onto  the  desired 
station  using  the  interpolants  specified  by  mdintr  .  The  and  derivatives  of 
the  basic  state  are  then  obtained  from  the  corresponding  derivatives  of  the 
interpolants. 

It  is  important  to  note  that  subroutine  bstate  uses  the  same  grid  points  chy(j)  at 
all  the  stations  specified  by  mdintr  even  when  pointd .  ins  maJces  the  disturbance 
grid  point  distribution  vary  from  station  to  station.  This  makes  it  easier  to  numer¬ 
ically  evaluate  the  derivatives  of  the  basic  state  with  respect  to  and  holding 
constant .  The  vector  chy(j)  may,  however,  be  different  the  next  time  that  bstate 
is  called  but  again  will  be  locally  held  constant  during  the  interpolation. 

Subroutine  bstate  reads  the  data  it  needs  for  interpolation  from  a  direct  access  file. 
This  data  access  method  enables  bstate  to  go  directly  to  the  appropriate  place  in  the 
file  and  read  in  the  desired  point  of  any  station  profile  without  having  to  first  read 
in  all  of  the  points  in  the  preceding  profiles.  The  direct  access  method  also  enables 
bstate  to  read  the  profiles  in  any  order  that  the  user  specifies,  even  from  last  to 
first.  Hence,  the  stability  investigation  can  begin  at  some  downstream  location  where 
it  is  often  easier  to  find  unstable  eigenvalues.  From  there,  the  analysis  can  proceed 
upstream  by  continuation  into  the  stable  region. 

The  FORTRAN  unformatted  direct  access  file  is  set  up  by  subroutine  gridbs  through 
calls  to  a  number  of  additional  routines^.  These  additional  routines  are:  bsname, 
bsiois,  bdiois,  bpstda,  and  scale. 

Two  of  these  routines  are  independent  of  the  basic  state  input  file  formats  and  are 
stored  in  file  util.f .  They  are:  bsname  and  bdiois.  util.f  also  contains  subroutine 
mdintr  because  mdintr  should  not  need  to  be  modified  for  different  flow  solver  inter¬ 
faces.  None  of  the  routines  in  util.f  use  any  variables  in  common. 

Subroutines  bsiois,  bpstda,  and  scale  are  stored  together  in  file  auz./ and  may  need 

^The  direct  access  file  would  not  really  have  to  be  created  if  each  input  file  itself  were  direct 
access,  but  this  is  rarely  the  case. 
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to  be  changed  when  developing  a  new  interface  to  a  different  basic  state  flow  solver. 
We  have  tried  to  circumvent  the  need  for  extensive  revisions  whenever  possible  by 
developing  an  interface  to  the  popular  PL0T3D  format.  Nevertheless,  not  every  flow 
solver  creates  PL0T3D  files. 

The  additional  routines  that  bstate  and  gridbs  use  to  help  with  the  data  interpo¬ 
lation  and  file  management  are  described  in  the  following  subsections. 


C.l  mdintr 

As  mentioned  above,  subroutine  mdintr  (MultiDimensionad  INTeRpolation)  deter¬ 
mines  which  stations  to  use  for  interpolating  the  basic  state  solution  in  the  variables 
and  It  also  computes  the  Lagrange  interpolants  and  their  derivatives  with 
respect  to  and 

The  Lagrange  interpolants  are  products  of  polynomicds  of  orders  ngk  and  ngj  in 
and  respectively.  Tentative  values  of  ngk  and  ngj  are  passed  as  arguments 
to  mdintr  .  However,  the  input  ngk  and  ngj  may  be  changed  if  they  exceed  the 
corresponding  basic  state  grid  dimensions  or  the  bounds  of  the  memory  allocated  for 
the  interpolants. 

Either  one  or  two-dimensional  interpolation  may  be  required  depending  upon  the 
dimensions  of  the  basic  state  solution. 

One-dimensional  interpolation  with  Lagrange  polynomials  will  be  used  in  the  direction 
for  which  the  grid  dimension  is  greater  than  one.  The  interpolant  in  the  other  direction 
(with  grid  dimension  1)  will  be  set  to  1.  Only  a  one-dimensional  search  scheme  is 
needed  in  this  case. 

Two-dimension<il  interpolation  with  Lagrange  polynomials  will  be  used  whenever  both 
grid  dimensions  are  greater  than  one.  In  this  case,  a  two-dimensional  search  scheme 
will  be  used  to  find  the  appropriate  place  in  the  basic  state  grid.  The  scheme  imple¬ 
mented  here  uses  a  variant  of  Newton-Raphson  iteration  to  solve  the  equations: 

=  i\RK,BJ)  (C.7) 

=  (HRK.RJ)  (C.8) 

for  the  basic  state  grid  indices  {RK,RJ)  corresponding  to  the  point  This 

scheme  may  not  always  work,  even  when  the  Jacobian  is  never  zero.  However,  it  is 
fast  and  in  practice  has  been  reasonably  robust. 

Once  the  real  values  RK  and  RJ  have  been  determined,  neighboring  stations 
(^^,^1),  k=l,2,, .  .,ngk,  j  =  l,2,...,ngj  are  selected  for  the  interpolation.  Tenta¬ 

tively,  nbk=ngk/2  points  below  and  ngk-nbk  points  above  ^p  in  the  basic  state  grid 
are  used  for  interpolating  in  : 

(it-1))  (C.9) 
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where  J{RK)  is  the  integer  obtained  by  eliminating  the  £raction«il  part  of  RK.  Sim¬ 
ilarly,  nbj»ngj/2  below  and  ngj-nbj  points  above  are  used  for  interpolating  in 

=  eiint{RJ)-nbj  +  {j-l))  (C.IO) 

The  stndat  addresses  of  all  these  stations  are  stored  in  the  integer  vector 
ida(l:nsrf)  .  The  total  number  of  points,  nsrf  *  ngk*ngj,  is  later  output  by 
subroutine  mdintr  . 

After  determining  which  stations  to  use  for  the  interpolation,  subroutine  mdintr  com¬ 
putes  the  interpolating  polynomials.  These  interpolants  are  taken  to  be  polynomials 
in  instead  of  {K,J)  because  the  basic  state  flow  usually  varies  smoothly  with 

while  it  may  not  with  {K,J).  Basic  state  grid  spacing  discontinuities  which 
are  often  present^  cause  discontinuities  in  the  derivatives  of  the  flow  with  respect  to 
{K,J).  Such  discontinuities  make  it  more  difficult  to  use  higher  order  interpolants. 
To  ensure  that  the  interpolants  in  2ure  well  defined^,  the  transformation  from 

the  basic  state  coordinates  {K,J)  to  the  disturbance  coordinates  is  restricted 

to  be  of  the  form: 


=  ^^(A'),A:  =  l,2,...,kdim  (C.ll) 

e  =  f(J),J  =  l,2,...,jdim  (C.12) 

where  kdim  and  jdim  are  the  dimensions  of  the  basic  state  grid  in  the  K  and  J 
directions,  respectively.  This  restriction  can  be  eliminated  by  interpolating,  say,  in 
the  arc  length  along  the  basic  state  grid  lines.  However,  the  restriction  has  not  yet 
been  removed. 

The  use  of  the  rectilinear  grid  in  equations  (C.ll)  and  (C.12)  above,  enables 

the  multi-dimensional  interpolation  to  be  done  by  composing  two  one-dimensionad 
interpolants; 


ngj 


F{e,e) 

(C.13) 

j  — a 

ngk 

fe=l 

(C.14) 

or 

F((\e)  = 

ngj  ngk 

j=l  fc=l 

(C.lo) 

^especially  in  the  files  which  the  AFWAL  PNS  code  creates  because  of  the  seemingly  arbitrary 
way  in  which  it  prints  out  data 

^See,  for  instance,  Lancaster  Ic  Salkauskas  (1988),  pp. 135-137. 
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where  F  is  an  arbitrary  fxmction  and 


Me)  = 
Me)  = 


(C.16) 

(C.17) 


The  values  of  tl^kie)  evaluated  at  4p  are  stored  in  bik(0,k)  .  The  ndk  derivative  of 
i^kie)  with  respect  to  evaluated  at  is  stored  in  bik(ndk,k)  .  Similarly,  the 
ndk  derivative  of  Me)  with  respect  to  ^  ,  evaluated  at  is  stored  in  bij  (ndj  ,  j)  . 
The  order  ndmx  of  the  highest  derivative,  mixed  or  otherwise,  is  passed  as  an  input  to 
mdintr  .  The  derivative  of  the  tensor  product  interpolant  with  respect 

to  ndk  times  and  ndj  times  is  stored  in  the  matrix  bi(idad,ia)  ,  where  ia  = 
ngj*(k-l)+j  and  idad  =  ndj*(2*ndmx-ndj+3)/2+ndkfor  0  <  ndj  <  ndmx  and  0  < 
ndk  <  ndmx-ndj 

All  of  the  station  data  is  interpolated  and  numerically  differentiated  using  the  matrix 
bi  .  If  the  body  coordinates  are  stored  in  stndat,  this  computes  the  body  slopes, 
body  curvature,  etc.  Control  is  then  returned  to  the  calling  program. 


C.2  bsname 

Subroutine  bsname  determines  the  names  bsfn(n)  (n=l,2,. .  ..nbsf)  of  the  nbsf 
basic  state  input  files.  It  also  creates  a  corresponding  direct  access  file  name  for  each 
basic  state  input  file.  The  other  subroutines  which  read  end  write  the  beisic  state 
information  may  then  use  as  many  direct  access  files  as  there  atre  baisic  state  input 
files.  However,  the  present  scheme  stores  all  of  the  information  in  one  direct  access 
file  only:  the  first  one. 

The  algorithm  that  bsname  uses  to  determine  the  name  bsfn(n)  of  input  file  n  is  as 
follows.  First,  it  checks  to  see  if  the  file  fort  .un  exists.  Here,  =  ubsf  (n)  is  the 
(integer)  unit  that  file  number  n  will  be  connected  to  for  input.  If  fort.un  exists, 
then  bsnauae  assigns  bsfn(n)  =  fort.un  .  The  user  can  enable  this  by  simply  soft- 
linking  the  input  files  to  the  corresponding  files  fort.un  .  If  fort.un  does  not  exist, 
then  bsname  prompts  the  user  with  a  message  to  input  the  corresponding  file  name. 
The  message  is:  "Enter  the  name  of  the  basic  state  ",  messag(lb:le) , " 
file  >  ",  where  messag(lb:le)  is  the  n‘*  substring  ending  with  a  semicolon  de¬ 
limiter  in  the  character  variable  messag  passed  to  bsnaune  by  gridbs^  . 

Once  bsname  has  determined  bsfn(n) ,  it  assigns  a  corresponding  name  to  daf  n(n)  us¬ 
ing  the  following  scheme.  If  bsfn(n)  =  fort  .ubsf  (n),  then  bsname  assigns  dafn(n) 
=  fort.uda^  .  If  bsfn(n)  ^  fort  .ubsf  (n),  then  bsname  strips  off  any  extension 

if  messag  =  ’grid;solution;par2uneters;’,  then  the  second  message  is  “Enter  the  name  of  the 
basic  state  solution  file  > 

^The  file  extension  ”uda”  stamds  for  “unformatted  direct  access” . 
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following  in  the  file  name  bsfn(n)  and  replax:es  it  with  "uda"  .  If  bsfn(ii)  does 
not  have  an  extension,  then  bsiiame  adds  the  extension  ”  .  uda”  .  In  any  of  these  cases, 
the  extension  for  each  file  might  not  be  unique.  To  avoid  this,  the  base  name  o/bsfn 
must  be  unique  in  order  for  the  corresponding  name  daifn  to  be  unique  .  This  is  not 
really  important  in  the  present  version  of  the  stability  codes  because  they  only  use 
one  direct  access  file  anyway. 

The  arguments  of  subroutine  bsname  axe:  nbsf,  ubsf  ,uda,bsfn,dafn,  and 
messag  .  ubsf  ^lnd  uda  are  integer  vectors  of  length  nbsf,  while  bsfn  and  daifn 
Me  chMacter  vectors  of  length  nbsf  .  messaig  is  a  chMacter  vMiable  with  nbsf  sub¬ 
strings,  each  substring  ending  in  a  semicolon  delimiter,  bsnauae  returns  to  the  calling 
routine  after  determining  all  of  the  names  of  the  files.  It  does  not  itself  open  or  close 
any  files. 

C.3  bsiois 

Subroutine  bsiois  reads  into  core  memory  cdl  the  basic  state  data  which  depend 
only  on  amd  Examples  of  such  data  Me:  number  of  points  in  the  profile; 
body  radius;  shock  stand-off  distance,  etc.  All  of  this  data  is  stored  in  the  Mray 
stndat (0:iistpar,oumsta),  where  nstpaor  is  the  number  of  station  parameters  and 
mxnsta  is  the  maximum  number  of  stations.  To  utilize  this  information,  the  interpo¬ 
lation  routines  all  require  that  it  be  recorded  in  the  following  elements  of  stndat: 

stndat  (0  ,k)  =  the  cumulative  number  of  points  in  all  the  profiles 
preceding  station  k,  plus  one.  The  user  must  have  bsiois 
compute  this  number.  Subroutine  bdiois  will  later  add  to 
stndat (0,k)  and  replace  it  with  the  direct  access  file  record 
number  for  the  beginning  of  profile  k  after  determining  how 
much  spaw:e  is  required  for  adl  of  the  problem  specific  data  and 
station  data. 

stndat  (1  ,k)  =  the  value  of  for  station  k  .  Must  be  set  by  sub¬ 
routine  bsiois  after  being  read  or  computed  from  the  basic 
state  input  files. 

stndat  (2  ,k)  =  the  number  of  grid  points  in  the  profiles  at  sta¬ 
tion  number  k  .  Must  be  set  by  subroutine  bsiois  after  being 
read  from  the  basic  state  input  files. 

stndat (3, k)  =  the  value  of  for  station  k  .  Must  be  set  by  sub¬ 
routine  bsiois  after  being  read  or  computed  from  the  baisic 
state  input  files. 
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Additional  information  such  as  the  Cartesian  coordinates,  body  radius,  or  shock 
stand-off  distance  may  also  be  stored  in  stndat  to  be  used  by  metric  or  btraois  . 
The  assignment  of  the  disturbance  coordinates  to  the  basic  state  grid  points  is  vital 
because  it  enables  the  interpolation  routines  to  find  their  way  around  in  the  basic  state 
grid.  This  may  require  that  bsiois  transform  the  basic  state  coordinates  (A',  J)  to 
the  disturbance  coordinates  if  they  are  not  the  same.  The  transformation  of  I 

to  must  be  made  by  subroutine  bpstda,  described  below,  as  it  reads  in  the  profiles 
and  writes  them  out  to  the  direct  access  file. 

Subroutine  bsiois  also  assigns  problem  specific  data  to  one-dimensional  vectors. 
There  is  one  vector  available  for  each  data  type  expected:  logicv(0  :nlp)  for  logical, 
charv(0:ncp)  for  character,  intv(-6:nip)  for  integer,  and  realv(0:nrp)  for  real. 
The  logical  vector  may  contain  flags  which  cure  set  by  the  baisic  state  flow  solver.  The 
character  vector  typically  contains  some  lines  of  a  title.  The  integer  vector  contains 
the  basic  state  grid  dimensions.  Finally,  the  real  vector  usually  contains  the  reference 
temperature  and  free-stream  Mach  number,  etc.  In  order  to  access  this  data,  the 
other  routines  require  that  bsiois  make  the  following  assignments. 


integer 

intv(-6)  = 
intv(-5)  = 
intv (-4)  = 
intv (-3)  = 
intv(-2)  = 
intv(-l)  = 
intv(  0)  = 
intv(  1)  = 
intv(  2)  = 
intv(  3)  = 


nip,  number  of  integer  variables, 
nip,  number  of  logical  variables, 
ncp,  number  of  character  variables, 
length  of  character  variables, 
nrp,  number  of  real  variables, 
number  of  station  parameters, 
number  of  basic  state  variables. 

K  dimension  of  basic  state  grid, 
total  number  of  stations  in  file. 

J  dimension  of  basic  state  grid. 


real 


scale  of  basic  state  grid. 

scale  of  basic  state  grid. 

scale  of  bcisic  state  grid. 

=  T  =  t  scale  of  basic  state  grid. 

realv(4+k)  =  scale  of  bcisic  state  solution  variable  number  k, 
k=l,2,. .  ..nbsvar,  where  nbsvar  =  intv(O)  is  the  number 
of  basic  state  variables. 


realv(l)  = 
realv(2)  = 
realv(3)  = 
realv(4)  = 
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C.4  scale 


After  all  of  the  station  data  has  been  read  in  from  the  original  basic  state  in¬ 
put  files,  subroutine  gridbs  calls  scale  .  Subroutine  scale  must  assigns  the 
scales,  scaldg(l:4)  ajid  scaldv(l:5),  respectively,  for  the  disturbance  coordinates 
and  default  disturbance  variables  (p,  T,  Here,  it‘  is  the 

velocity  component  in  inertial  Cartesian  coordinates,  scale  must  also  set  the  ratios 
scalbgCl  :4)  of:  the  disturbance  coordinate  scales  used  by  the  basic  state  solution  to: 
the  corresponding  scales  assigned  for  the  stability  analysis.  Finally,  it  also  must  set 
the  ratios  scalbv(l  :nbsvar)  of;  the  basic  state  input  variables  to:  the  correspond¬ 
ing  scales  assigned  for  the  stability  analysis.  The  basic  state  scales  set  in  subroutine 
bsiois  and  stored  in  the  direct  access  file  are  usually  used  here^.  The  ratios  are  later 
used  by  subroutines  mdintr  and  bstate  to  properly  rescale  the  basic  state  input  for 
the  stability  analysis. 


C.5  bdiois 

The  first  time  that  subroutine  gridbs  calls  bdiois,  bdiois  copies  all  of  the  sta¬ 
tion  data  and  problem  specific  data  to  the  top  of  the  direct  access  file.  Subroutine 
gridbs  calls  bdiois  prior  to  subroutine  bpstda  (see  below)  because  bdiois  com¬ 
putes  stndat(0,k) ,  k=l,2,. . intv  (2),  the  direct  access  file  record  number  for  the 
beginning  of  each  basic  state  profile  ^ 

The  next  time  that  the  program  is  run,  the  direct  access  file  will  already  exist.  In 
this  case,  subroutine  gridbs  will  still  call  bsname  to  get  the  name(s)  of  the  direct 
access  file(s).  However,  gridbs  need  not  call  subroutines  bsiois  and  bpstda  to  get 
the  stations  data  and  profiles.  Instead,  subroutine  gridbs  will  call  bdiois  to  read 
the  station  data  from  the  direct  access  file.  Subroutine  bstate  will,  as  before,  read 
the  profiles  from  the  direct  access  file. 


C.6  bpstda 

Subroutine  bpstda  copies  all  of  the  basic  state  input  profiles  to  the  direct  access  file. 
All  of  this  data  must  be  scaled  as  specified  by  the  assignment  statements  to  realv  in 
subroutine  bsiois  .  The  profiles  must  be  written  to  the  unformatted  direct  access 
file  using  statements  similar  to  the  following: 

do  20  k=l,intv(2) 

*When  a  particular  variable  really  has  no  scale,  its  scale  can  be  set  arbitrarily  as  long  as  the 
scale  is  consistent  with  the  dimensionless  parameters  . 

^The  record  number  is  computed  and  stored  before  the  array  stndat  is  copied  to  the  direct  access 
file. 
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do  10  ii=l,nint(stndat(2,k)) 


C 

C  Read  in  the  basic  state  variables  qb(nbsv)  at  point  n  in  profile  k. 

C  Compute  at  point  n  in  profile  k. 

C 

writeCudad)  .rec^key+n-l)  ^^,(qb(m)  ,m=l,nbsv) 

10  continue 

20  continue 

where  intv(2)  is  the  number  of  basic  state  profiles,  stndat(2,k)  is  the  number  of 
points  in  profile  k,  and  nbsv  =  intv(O)  is  the  number  of  basic  state  variables.  Again, 
the  integer  vector  intv  must  be  previously  set  in  subroutine  bsiois  . 

Once  the  direct  access  file  dafn(l)  has  been  created,  the  input  bcisic  state  files  may 
be  closed.  The  only  routines  that  will  subsequently  execute  will  be  the  interpolation 
routines  mdintr  and  bstate  .  They  will  only  use  the  information  in  the  direct  access 
file. 
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